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ABSTRACT 


Loss  of  constraint,  which  has  been  obser\'ed  in  various  laboratory 
specimens,  is  caused  by  the  interaction  of  crack  tip  plaisticity  with  free 
boimdaries.  When  this  occurs,  toughness  is  geometry  dependent,  and 
classical  single-parameter  fracture  mechanics  is  no  longer  valid.  A  number 
of  researchers  have  addressed  constraint  loss  by  introducing  a  second 
parameter  into  fracture  mechanics  methodologies.  There  are  instances, 
however,  when  even  two-parameter  fracture  tiieory  breaks  down. 

This  report  addresses  the  limitations  of  two-parameter  fracture  me¬ 
chanics.  We  performed  an  asymptotic  analysis  of  the  general  power  series 
representation  of  the  crack  tip  stress  potential  in  an  elastic  plastic  material 
that  obeys  a  Ramberg-Osgood  constitutive  law.  Expansion  of  the  power 
series  over  a  substantial  niomber  of  terms  yields  orily  tiiree  independent 
coefficients  for  low-  and  medivim-hardening  materiak.  The  first 
independent  coefficient  is  uniquely  related  to  the  J-integral,  as  Hutchinson, 
Rice  and  Rosengren  discovered  over  25  years  ago.  The  second  and  diird 
independent  coefficients,  K2  and  K4,,  are  a  function  of  geometry  and 
loading  level.  A  two-parameter  theory  implies  that  the  crack  tip  stress 
fields  have  two  degrees  of  freedom,  but  the  asymptotic  analysis  implies 
that  three  parameters  are  required  to  characterize  near-tip  conditions.  Thus 
two-parameter  fracture  theory  is  a  valid  engineering  model  only  when 
there  is  an  approximately  imique  relationship  between  K2  and  K4. 

We  performed  elastic-plastic  finite  element  analyses  on  several  geome¬ 
tries  and  evaluated  K2  and  X4  as  a  function  of  deformation  level.  A  refer¬ 
ence  two-parameter  solution  (which  gives  a  unique  relationship  between 
K2  and  X4)  was  pro\dded  by  the  modified  boundary  layer  (MBL)  geometry 
with  an  applied  Mode  I  stress  intensity  factor,  Ki,  and  various  T  stress  vzJ- 
ues.  The  MBL  configuration  simulates  a  crack  in  an  infinite  body.  The  fi¬ 
nite  geometries  we  considered  w’ere  the  single  edge  notched  bend  (SENB), 
single  edge  notched  tension  (SENT),  double  edge  notched  tension  (DENT), 
and  center  cracked  tension  (CCT)  specimens.  Crack  length-to-width  (a/W) 
ratios  of  0.1,  0.5  and  0.9  were  studied. 

Resiilts  indicate  that  the  near  tip  stresses  in  all  but  the  deeply  cracked 
SENB  (a/W  =  0.5,  0.9)  and  SENT  (a/W  =  0.9)  lend  fiiemselves  to  a  two-pa¬ 
rameter  characterization.  However,  the  deeply  cracked  SENB  and  SENT 
specimens  maintain  a  high  level  of  constraint  to  relatively  large  deforma¬ 
tion  levels.  Thus  single-parameter  fracture  mechanics  is  fairly  robust  for 
these  high  constraint  geometries,  but  two-parameter  theory  is  of  little 
value  when  constraint  loss  eventually  occurs.  On  the  other  hand,  in  low 
constraint  geometries  such  as  the  CCT  panel  and  shallow  notched  SENB 
specimens,  the  deparfiire  from  single-parameter  theory  is  almost  immedi¬ 
ate,  but  the  deformation  evolves  according  to  two-parameter  theory  (i.e., 

K2  and  K4  maintain  a  nearly  unique  relationship). 
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CHAPTER  I 


INTRODUCTION 


Since  the  introduction  of  what  is  known  today  as  single  parameter 
linear  elastic  fracture  mechanics  (LEEM)  in  the  early  to  mid  twentieth  century 
[1-6],  the  focus  of  research  in  this  discipline  has  had  two  primary  objectives. 
First,  it  was  necessary  to  apply  the  methodology  known  as  fracture  mechanics 
to  real  world  problems.  Standards  have  evolved  to  ensure  that  this  was  done 
in  a  systematic  way.  Second,  within  the  original  methodology,  simplifying 
assumptions  (e.g.  linear  elastic,  single  parameter,  boundaries  far  removed 
from  crack-tip,  etc.)  were  made  to  make  the  problems  mathematically  more 
tractable.  As  issues  such  as  plasticity,  crack  growth,  and  finite  geometries 
came  to  being,  the  theory  behind  the  original  methodology  had  to  be  changed 
to  include  effects  that  result  as  a  consequence  of  these  issues. 

Through  the  nineteen-sixties  and  early  seventies,  much  attention  was 
focused  on  the  aspect  of  crack-tip  plasticity.  Several  works  [7-9]  were 
completed  which  extended  the  idea  of  single  parameter  fracture  mechanics  to 
materials  that  exhibit  nonlinear  constitutive  relatior\s. 

The  combination  of  improved  computer  hardware  and  computerized 
finite  element  techniques  in  the  seventies  to  mid-eighties  enabled  numerical 
modeling  of  elastic  and  limited  plastic  crack-tip  problems.  It  was  at  this  time 
that  the  previously  pure  analytic  techniques  could  be  compared  with 
numerical  results. 

The  increase  in  both  speed  of  central  processing  units  (CPU's)  and 
accessible  memory  of  computers  in  the  late  eighties  and  early  nineties  made  it 
possible  to  numerically  model  complex  geometries  and  include  plasticity 
aspects.  Although  the  theoretical  methodology  seemed  to  benefit  from  the 
increased  computational  capabilities  developed  through  the  eighties,  the 
majority  of  the  test  standards  are  based  on  the  original  analytical  work  of  the 
forties  through  the  sixties  [10-12]. 

Constraint  has  become  an  important  issue  in  the  transfer  of  fracture 
The  Journal  Model  is  Theoretical  and  Applied  Fracture  Meclmnics. 
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mechanic  laboratory  results  to  structural  applications.  The  loss  of  constraint 
resulting  from  the  development  of  plastic  fields  and  their  interaction  with  near 
tip  finite  boundaries  found  in  many  laboratory  specimens  and  structures  has 
led  to  studies  whose  objectives  were  to  quantify  the  loss. 

The  original  theoretical  work  in  fracture  mechanics  avoided  these 
issues  by  assuming  a  field  characterized  by  a  single  parameter  in  which 
structural  boundaries  were  far  removed  from  the  crack-tip.  Most  recently, 
attempts  have  been  made  to  accoimt  for  the  constraint  issue  by  introducing  a 
second  parameter  to  the  methodology  [13-20].  Several  studies  [21-30]  had  for 
their  objectives,  the  idea  of  quantifying  the  limits  of  the  single  parameter 
methodology,  but  to  the  author's  knowledge,  no  attempt  has  been  made  to 
determine  the  limits  of  the  two-parameter  methodology.  It  is  the  author's 
intent  to  quantify  the  applicability  of  the  present  two-parameter 
methodologies  by  means  of  higher  order  terms  in  the  general  power  series 
representation  of  the  stress  potential. 

1.1  Present  Status  of  Problem 

The  original  work  of  Williams  [4]  was  based  on  an  asymptotic  analysis 
of  the  general  power  series  representation  for  the  stress  potential  near  the  tip 
of  an  infinitely  sharp  notch, 

x[r,  0,  A)  =  sin|(A  -1-1)0} +0^  cos|(A  +1)0} 

+ sin|(A  - 1)0} + cos{(A  -  l)0}j 

where  “  indicates  the  dimensional  quantity.  This,  combined  with  a  linear 
constitutive  relation,  was  used  with  small  strain  theory  to  obtain  the  pseudo 
eigenvalue  A.  The  boundary  conditions  for  the  crack  problem  are  satisfied  for 
the  infinite  sequence  of  values  A  =  n/2  where  n  =  1, 2, 3, ....  Williams  was  able 
to  determine  the  series  representation  of  the  stress  field  from  the  potential  as 
given  by 
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(1.2) 


where  G^.  are  angular  functions  and  are  related  to  the  multiplicative 

constants  Ci.  Williams  was  able  to  characterize  the  stress  field  up  the 
multiplicative  constants 

At  the  same  time  Williams  was  performing  his  asymptotic  analysis, 
Irwin  [5]  was  also  analyzing  crack-tip  fields  in  linear  media.  Irwin's  approach 
utilized  the  complex  stress  potential  of  Westergaard  [6], 

F  =  ReZ(z)+yImZ(z).  (1.3) 


By  taking 


ZU)= 


(1.4) 


the  stress  field  for  a  central  straight  crack  located  on  the  x-axis  in  an  infinite 
plate  with  a  biaxial  field  of  tension,  o’,  could  be  determined.  The 
corresponding  Cartesian  components  for  the  first  symmetric  term  in  the  near 
tip  stress  fields  for  this  and  Williams  method  are  given  as 


Although  the  angular  variation  of  the  crack-tip  stress  fields  found  by  Irwin 
were  identical  to  those  obtained  by  Williams,  Irwin  was  able  to  relate  the  first 
multiplicative  constant  in  the  stress  potential  to  the  far  field  boundary 
conditions.  This  constant  known  today  as  the  stress  intensity  factor  (JC)  is  the 
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basis  for  single  parameter  linear  elastic  fracture  mechanics  and  is  given  here  as 
K  =  a^Hta  for  the  central  crack  in  an  infinite  plate. 

Since  the  initiation  of  LEFM,  the  framework  for  a  multi-parameter 
methodology  was  available  but  practitioners  chose  to  ignore  the  higher  order 
terms  in  the-  stress  field.  Irwin  was  the  first  to  notice  empirically  the 
importance  of  the  higher  order  terms.  When  comparing  photo  elastic 
representations  of  stress  fields  to  his  analytic  predictions,  a  discrepancy  was 
foimd  in  the  stress  parallel  to  the  crack.  This  disaepancy,  as  shown  by  Sih 
[31]  theoretically,  is  now  known  as  the  T-stress.  It  is  the  second  term  in  the 
Williams  general  power  series  expansion.  This  term  gives  zero  contribution  to 
the  opening  mode  stress  but  exists  as  a  constant,  hence  no  angular  variation  in 
the  stress  parallel  to  the  crack  plane. 

Hutchinson  [7]  used  the  asymptotic  approach  of  Williams  in  his 
analysis  of  a  Ramberg-Osgood  power-law  hardening  material. 


Here  v  is  poisson's  ratio,  is  the  yield  stress,  Cg  is  the  effective  stress,  S  is 
the  deviatoric  stress,  Bq  is  the  2%  offset  yield  strain,  and  a  is  a  fitting 
parameter.  Hutchinson  was  able  to  numerically  determine  the  angular 
variation  of  the  first  term  in  the  series  representation  of  the  stress  field  as  it  is 
derived  from  the  stress  potential  and  given  by 

a-j(r,e)  =  A,o;J(e)/'.  (1.7) 

Concurrently,  Rice  and  Rosengren  [8]  were  using  the  asymptotic  approach  in 
power-law  hardening  materials.  The  first  multiplicative  constant  was  related 
to  the  far  field  boundary  conditioris  through  Rice's  [9]  path  independent  /- 
integral  and  is  given  by 


A^  =  ao 


J 


j 


(1.8) 
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This  work  extended  the  methodology  of  single  parameter  frachire  mechanics 
from  linear  elastic  materials  to  those  that  behave  according  to  deformation 
type  plasticity. 

The  main  assumption  of  this  theory  was  the  description  of  the  crack-tip 
stress  field  by  a  single  parameter,  namely  the  first  multiplicative  factor  in  the 
general  power  series  representation  of  the  stress  potential. 

The  first  attempt  to  accoimt  for  higher  order  asymptotic  terms  in 
power-law  hardening  materials  was  provided  by  Yaochen  and  Wang  [13]. 
They  carried  out  a  complete  two  term  asymptotic  analysis  and  were  able  to 
determine  using  a  fourth  order  numerical  Runge-Kutta  method  the  angular 
variation  of  the  second  term  in  the  stress  potential.  Yaochen  et  al.  were  able  to 
determine  the  second  multiplicitive  constant  by  comparing  the  two  term 
analytic  solution  to  finite  element  results  obtained  by  Shih  and  German  [24] 
and  Needleman  and  Tvergaard  [25]  for  a  single  edge  notched  bend  bar 
(SENB),  single  edge  crack  panel  (SECP),  and  a  center  cracked  panel  (CCP) 
having  a/ W  ratios  fi-om  0.5  to  0.9  and  hardening  exponents  of  3, 10,  and  «>. 

Yaochen  et  al.  concluded  that  the  second  coefficient  of  the  asymptotic 
field  changes  with  specimen  geometry  and  development  of  the  plastic  zone. 
As  a  direct  result,  the  second  multiplicative  factor  Kz  was  proposed  as  a 
second  parameter,  necessary  in  some  instances  to  characterize  the  stress  field 
at  the  crack-tip.  In  so  doing,  Yaochen  et  al.  introduced  the  idea  of  two- 
parameter  fracture  mechanics. 

Sharma  and  Aravas  [14]  expanded  on  the  work  of  Yaochen  et  al.  by 
solving  for  the  second  asymptotic  term  in  both  plane  strain  and  plane  stress. 
Sharma  et  al.  solved  the  problem  in  a  different  manner  then  that  of  Yaochen. 
The  two  in-plane  equilibrium  equations  and  the  three  linear  strain 
displacement  relations  were  solved  using  finite  elements.  The  asymptotic 
results  were  verified  with  full  field  finite  element  calculations.  The  full  field 
calculations  were  based  on  the  boundary  layer  method  outlined  by 
McMeeking  [21],  and  Rice  and  Tracey  [32].  In  their  work,  they  applied 
displacements  corresponding  to  the  first  Williams  term  to  the  armular 
boimdary  of  the  crack-tip  region.  The  plastic  fields  that  result  correspond  to 
those  which  would  be  expected  at  the  tip  of  a  crack  under  high  constraint  for 
which  all  boundaries  were  far  removed.  These  results  have  become  known  as 
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the  small  scale  yielding  solution  (SSY).  The  plane  strain  results  obtained  for 
the  second  as5nnptotic  term  were  identical  to  those  of  Yaochen. 

O’Dowd  and  Shih  [15,16]  address  what  has  become  known  as  the 
difference  field  in  an  alternate  manner.  O'Dowd  et  al.  argue  that  the  crack-tip 
stress  field  could  be  expressed  as 

ay(r,e)=(f^{r,e)+QSij.  (1.9) 

Here  the  stress  field  has  been  normalized  by  They  show  the  difference 
field  represented  here  by  Q  to  be  hydrostatic  in  nature  in  the  forward  sector  of 

the  crack  out  to  a  normalized  radius  r  =  of  about  ten.  The 

magnitude  of  Q  is  determined  by  taking  the  difference  in  the  opening  mode 
stress  on  the  crack  plane  between  a  finite  body  and  that  for  the  SSY  solution. 
The  difference  between  these  fields  is  taken  at  a  predetermined  location  in  the 
forward  sector  generally,  0  =  0  and  r  =  2.  Although  strictly  speaking  this  is 
not  theoretically  rigorous,  it  seems  to  be  somewhat  robust. 

Xia,  Wang  and  Shih  [17]  and  Yang,  Chao  and  Sutton  [18,19]  expand  the 
asymptotic  analysis  for  a  power  hardening  material  to  the  fifth  term.  In  each 
study,  the  strain  hardening  exponent  was  varied  to  determine  which  term 
would  contain  the  first  elastic  field.  In  the  analysis,  only  three  amplitudes 
could  be  prescribed  independently.  For  the  appropriate  choice  of  the 
independent  amplitudes,  it  was  possible  to  reproduce  near  tip  fields 
representative  of  a  broad  range  of  crack-tip  constraints.  The  higher  order 
terms  obtained  in  the  studies  collectively  describe  an  approximately  spatially 
uniform  hydrostatic  stress  field.  These  results  lend  support  to  the  suggestion 
given  by  O'Dowd  and  Shih  [15,16]  that  Q  could  be  introduced  as  a  second 
fracture  parameter  that  quantifies  the  near  tip  triaxiality. 

Anderson  and  Dodds  [33],  Dodds,  Anderson  and  Kirk  [34],  and 
Anderson  and  Dodds  [35]  use  a  micro  mechanics  approach  to  address  the 
constraint  issue.  Unlike  the  afore  mentioned  work  generally  based  on 
continuum  mechanics,  Anderson  et  al.  base  their  model  on  a  local  criterion  for 
cleavage  failure.  The  model  assumes  that  cleavage  fracture  is  controlled  by 
the  volume  of  material  at  the  crack-tip  that  is  subject  to  a  high  stress  level;  the 


larger  the  stressed  volume,  the  higher  the  probability  that  cleavage  will 
initiate  from  a  critical  microstructural  feature.  The  loss  of  constraint  in  finite 
bodies  tends  to  reduce  the  volume  over  which  the  higher  stresses  preside.  As 
a  result,  one  finds  it  necessary  to  apply  a  higher  /  value  to  a  finite  body  to 
obtain  the  same  stressed  volume  foimd  in  an  infinite  body  at  a  lower  J  value. 

The  ratio  of  the  finite  body  /  value  to  the  corresponding  infinite  body 
value  defined  as  Jo  is  determined  by  the  ratio  of  volumes  contained  within 
similar  principle  stress  contours.  This  ratio  can  be  thought  of  as  a  scaling 
factor  that  enables  one  to  obtain  a  fracture  toughness  from  one  size  specimen 
and  then  scale  it  so  it  corresponds  to  that  which  would  be  obtained  from  an 
infinite  body  or  some  other  geometry.  This  approach  has  enabled  the 
correction  of  fracture  data  which  was  affected  by  the  loss  of  constraint 
common  in  many  laboratory  specimens. 

Why  should  one  be  interested  in  the  deformation  limits  on  two- 
parameter  fracture  mechanics?  The  stress  field  characterization  by  a  single 
parameter  breaks  down,  for  example,  when  free  boundaries  impinge  on  the 
crack-tip.  The  reduced  crack-tip  constraint  associated  with  these  free 
boundaries  results  in  a  lowering  of  the  crack-tip  stresses.  Data  presented  by 
Sorem  [36]  shows  an  increase  in  the  apparent  elastic-plastic  fracture  toughness 
(/c)  for  specimens  with  shallow  cracks  when  tested  in  the  trai^sition  region. 

In  the  transition  region,  the  fracture  mechanism  is  changing  from 
fracture  by  cleavage  which  is  predominately  stress  controlled  to  that  by 
ductile  rupture  which  tends  to  be  strain  controlled.  The  apparent  increase  in 
toughness  results  from  the  loss  of  crack-tip  constraint.  Cleavage  failure  occurs 
when  the  stress  over  some  small  but  measurable  crack-tip  volume  attains  a 
value  sufficient  to  trigger  the  failure  mechanism.  Therefore,  for  structures  that 
exhibit  lower  crack-tip  constraint,  the  applied  loading  must  be  greater  to 
produce  the  same  critical  crack-tip  stresses.  The  associated  increase  in 
toughness  is  analogous  to  the  increase  observed  when  changing  from  plane 
strain  to  plane  stress. 

On  the  opposite  end  of  the  transition  region  and  on  the  upper  shelf 
were  ductile  rupture  governs  failure,  Joyce  [37]  has  shown  that  the  loss  of 
constraint  associated  with  near  tip  finite  boundaries  has  little  influence  on  the 
toughness  measured  at  the  initiation  of  ductile  crack  growth  (Jic).  The  data 
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does  indicate  that  the  loss  of  constraint  influences  the  shape  of  the  R-curve,  or 
more  specifically  the  apparent  tearing  modulus.  The  data  shows  a  somewhat 
linear  relation  between  the  tearing  modulus  (measured  at  a  fixed  crack 
extension)  and  constraint  measured  in  terms  of  the  Q  parameter  introduced  by 
O’Dowd  and  Shih  [15,16], 

The  initiation  of  ductile  fracture  occurs  as  a  result  of  void  nucleation 
around  second  phase  particles  and  is  followed  by  void  growth  until  the  voids 
coalesce.  Argon,  Im  and  Safoglu  [38]  have  argued  that  the  initiation  of  voids 
around  second  phase  particles  occur  at  a  critical  combination  of  the 
hydrostatic  and  effective  stress.  Since  the  loss  of  constraint  has  little  effect  on 
the  value  of  fic,  one  could  postulate  the  effective  stress  which  depends  on  the 
maximum  differences  in  the  principle  stresses  is  more  important  than  the 
hydrostatic  stress  in  the  void  nucleation  process,  A  simple  examination  of 
Mohr’s  circle  demonstrates  that  the  effective  stress  can  remain  constant  for  any 
value  of  the  hydrostatic  component  of  the  stress  field. 

Rice  and  Tracy  [39]  have  also  shown  that  void  growth  is  dependent  on 
the  principle  stresses  and  strain  rate.  Here  the  absolute  magnitude  of  the 
principle  stress  components  are  decreased  with  a  reduction  of  the  hydrostatic 
stress.  Since  the  void  growth  process  appears  to  be  dependent  on  the 
magnitude  of  the  principle  stresses,  this  could  account  for  the  dependence  of 
the  R-curve  shape  on  the  degree  of  crack  tip  constraint  associated  with  a  given 
specimen  geometry.  It  would  appear  that  the  effective  stress  ahead  of  the 
crack-tip  is  sufficient  (even  with  small  decreases  in  the  hydrostatic  stress 
associated  with  loss  of  constraint)  to  initiate  void  nucleation  and  therefore 
void  growth  and  coalescence  are  the  important  stages  of  ductile  fracture  as  it 
relates  to  crack-tip  constraint, 

1.2  Research  Objective 

As  is  apparent,  fracture  toughness  measures  are  dependent  on  the 
degree  of  constraint  to  which  a  crack-tip  is  subjected.  The  current  method  for 
dealing  with  the  issue  of  constraint  is  to  introduce  a  second  parameter  which 
corrects  the  stress  level  on  the  crack  plane  to  agree  with  that  measured 
numerically.  The  second  parameter,  along  with  the  first  represent,  a  two 
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dimensional  failure  locus.  Stress  states  that  correspond  to  points  within  the 
locus  are  safe  while  points  outside  the  locus  are  subject  to  fracture. 

There  are  three  two-parameter  methodologies  at  present.  Each  of  the 
methodologies  use  the  /-integral  to  quantify  the  deformation  of  fracture.  The 
difference  in'  the  three  methodologies  is  in  the  definition  of  the  second 
parameter.  The  second  parameter  in  one  sense  or  another  is  a  measure  of  the 
hydrostatic  shift  in  the  crack  plane  stresses  associated  with  changes  in  the 
degree  of  constraint  relative  to  that  of  high  constraint  or  the  small  scale 
yielding  model.  One  approach  uses  the  value  of  the  elastic  T-stress  just 
outside  the  plastic  zone  to  measure  constraint,  while  another  method  uses  the 
difference  between  the  normal  stress  to  the  crack  plane  for  small  scale  yielding 
and  the  stress  measiued  in  a  finite  body  at  a  fixed  location  on  the  crack  plane. 
This  difference  is  known  as  the  Q -stress.  Although  the  Q-stress  is  not 
theoretically  rigorous  in  nature  it  appears  to  be  somewhat  robust  in  its 
application.  The  third  two-parameter  methodology  is  based  on  the  second 
coefficient  in  the  general  power  series  expansion  of  the  stress  field.  This 
method  although  being  more  theoretically  rigorous  suffers  because  it  is  only 
applicable  for  a  Ramberg-Osgood  constitutive  law.  The  definition  of  the 
second  parameter  is  in  itself  not  that  important  as  long  as  it  is  applied  in  a  self 
consistent  manner.  However,  the  T-stress  methodology  breaks  down  with 
increased  plasticity. 

The  objective  of  this  proposed  dissertation  research  is  to  evaluate  the 
deformation  limits  on  two-parameter  fracture  mechanics  as  defined  by  higher 
order  asymptotic  terms  in  nonlinear  crack-tip  stress  fields  that  are  influenced 
by  geometric  and/or  material  constraints. 

1.3  Proposed  Research  Program 

The  deformation  limits  of  two-parameter  elastic-plastic  fracture 
mechanics  will  be  studied  by  comparing  the  amplitudes  in  the  general  elastic- 
plastic  power  series  representations  for  stresses  within  the  plastic  zones  for 
two-parameter  plastic  fields  with  those  for  multi-parameter  plastic  fields. 

The  two-parameter  plastic  field  which  will  be  used  for  this  study  is 
based  on  the  solution  of  the  modified  boundary  layer  (MBL).  Consider  a 
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cracked  body  in  which  the  plastic  zone  is  well  contained  within  an  elastic 
region.  If  the  plastic  zone  is  small  relative  to  any  physical  dimensions,  the 
stresses  in  an  annulus  outside  the  plastic  zone  can  be  characterized  by  the 
general  elastic  power  series  representation  developed  by  Williams  [4].  The 
stresses  in  an  annulus  within  the  plastic  zone  but  outside  the  finite  strain  zone 
can  also  be  characterized  by  the  general  elastic-plastic  power  series 
representation  developed  by  Hutchinson  [7],  Rice  and  Rosengren  [8].  If  we 
assume  initially  that  the  conditions  in  the  far  field  lend  to  the  characterization 
of  the  stresses  in  the  elastic  annulus  by  the  first  two  terms  in  the  Williams 
solution,  then  the  fields  in  the  plastic  zone  by  definition  are  only  dependent 
on  the  same  two  parameters.  The  amplitudes  of  the  elastic-plastic  power 
series  themselves  are  dependent  on  these  two  parameters  also.  The 
amplitudes  in  the  elastic-plastic  series  will  be  determined  by  matching  the 
series  with  the  numerical  results  obtained  for  the  plastic  zone  in  the  MBL  for  a 
given  pair  of  amplitudes  in  the  elastic  power  series.  Rice  [9]  has  shown  the 
first  amplitude,  isTj,  in  the  elastic-plastic  series  is  completely  determined  using 
the  /-integral.  Yang,  Chao  and  Sutton  [19]  have  shown  that  there  are  only  two 
independent  amplitudes  (K2,K^)  remaining  in  the  first  five  terms  of  the  series. 
The  third  and  fifth  amplitudes  depend  on  the  independent  amplitudes  in  the 
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following  marmer  Therefore  a  unique  relation  between 
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the  two  independent  amplitudes  can  be  established  as  a  function 


of  deformation  level  since  the  plastic  fields  are  two-parameter  fields  and  the 
first  amplitude  is  determined  completely  by  the  /-integral. 

The  relationship  between  K2  and  for  the  two-parameter  plastic  field 
will  be  compared  with  a  similar  relationship  obtained  for  several  finite  bodies. 
The  characterization  of  the  stress  field  in  the  finite  bodies  can  depend  in 
general  on  a  number  of  parameters.  Comparison  of  the  K2  and  relation 
for  the  modified  boundary  layer  with  that  for  the  finite  bodies  as  shown 
schematically  in  Fig.  1.1  will  indicate  when  a  two-parameter  characterization 
of  the  stress  field  is  no  longer  appropriate.  It  is  possible  to  construct  a 
difference  field  between  the  two-parameter  results  and  the  finite  body  results 
as  given  by 
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where 


(1.11) 

When  the  error  defined  by  this  relation  is  greater  then  a  set  value  the  two- 
parameter  solution  is  deemed  no  longer  acceptable.  The  procediures  which 
will  be  used  to  construct  this  comparison  are  outlined  in  the  following  section. 
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CHAPTER  II 
PROCEDURES 


The  following  sections  contain  the  procedures  that  were  followed  to 
obtain  the  data  necessary  to  compare  the  two-parameter  fields  with  those  of 
the  finite  geometries.  The  constitutive  models  for  the  analytic  and  numerical 
analysis  are  presented.  The  procedures  used  to  solve  the  analytic  problem 
are  given  as  well  as  the  niunerical  methods  used  to  determine  the  finite 
geometry  results. 

2.1  Constitutive  Models 

The  total  strain  at  the  tip  of  a  crack  in  an  elastic-plastic  material  is 
generally  dominated  by  the  plastic  component.  When  unloading  is 
insignificant,  deformation  plasticity  has  proved  to  be  acceptable  in  aack-tip 
modeling.  In  the  present  work,  deformation  plasticity  is  used.  The 
constitutive  law  used  in  the  analytical  study  is  based  on  Ramberg-Osgood 
behavior.  The  material  for  numerical  studies  is  based  on  a  linear  plus  power- 
law  behavior  (Ramberg-Osgood).  The  transition  between  the  two  segments  is 
elliptic  in  nature.  In  the  analytic  study,  the  elastic  strains  in  the  Ramberg- 
Osgood  relation  do  not  enter  the  solution.  Hence,  the  dominate  strain  at  the 
crack-tip  is  assumed  to  be  power-law  in  nature.  Since  the  region  of  interest  is 
near  the  crack-tip,  both  the  Ramberg-Osgood  and  the  linear  plus  power-law 
give  equivalent  results.  The  use  of  the  linear  plus  power-law  material  is 
necessary  in  the  modified  boxmdary  layer  analysis  to  prevent  premature 
yielding  of  material  in  the  far  field  as  well  as  material  next  to  the  boundary. 

2.1.1  Analytical  Constitutive  Model 

The  Ramberg-Osgood  total  strain  is  assumed  to  be  the  sum  of  the 
elastic  and  plastic  components  as  given  by 
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pTotal  p  Bostic  pPlasiic 

_y _ __y _ _ 

^0  ^0  ^0 


(2.1) 


Unless  otherwise  stated,  ^  indicates  the  full  dimensional  quantity  while  the 
absence  of  the  "  indicates  a  normalized  value.  The  elastic  strain  is  assumed  to 
follow  Hooke's  law  and  is  given  as 


pEUutie 
_«/ _ 

£o 

The  plastic  strain  is  given  by  the  associative  flow  rule 

^Plastic  ,  „  -Plastic  S.. 
3/ _ _  J  0  ^ _ U_ 

£0  lata 

e  o  0 
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where  is  the  equivalent  effective  plastic  strain  given  by 


•^Plastic 
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^Plastic  Plastic 

2  c..  c.. 
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(2.4) 


and  dg  is  the  effective  stress  (Mises  stress)  determined  from 
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where  S^j  is  the  deviatoric  component  of  stress  written  as 
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Substitution  of  Eq.  (2.2-2.6)  into  Eq.  (2.1)  results  in 
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Noting  that  the  normalized  hydrostatic  stress  can  be  expressed  as 
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and  the  total  deviatoric  strain  is  given  by 


^Total  -  ^Total 

^ _ __y _ ilzn2L_5 

£o  e  3  £  y’ 


The  deviatoric  stress /strain  relation  can  be  expressed  as 


pToud 

ij  - 


r*  ^ 

-5  Cr  7 Plastic  S.. 


2  O’  e  O’ 
e  0  J  0 


(2.10) 


and  the  equivalent  effective  deviatoric  stress/strain  relation  is  obtained  by 
taking  the  inner  product  of  the  above  relation  to  give 


^Total  2  3<^  pPIastic  - 

£ - =  ±  (l  +  v)+-^^ - 

So  3  ^  ^  2a  £  CTo 

e  0  . 


where  the  total  effective  deviatoric  strain  is  given  by 


(2.11) 


-Total  .  e^otal  ^Total 
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(2.12) 


The  equivalent  plastic  strain  is  assumed  to  obey  a  power-law  relation 

given  by 
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(2.13) 


where  n  is  the  strain  hardening  exponent  and  a  is  a  fitting  parameter. 
Combining  the  above  expressions  results  in  the  following  expression  for  the 
total  deviatoric  Ramberg-Osgood  strain  in  terms  of  stresses 

J— =  (i  +  v)+-a -£■  (2.14) 


This  constitutive  relation  is  used  in  the  asymptotic  analysis. 


2.1.2  Numerical  Constitutive  Model 


The  constitutive  model  used  for  the  numerical  study  in  the  present 
research  is  a  linear  plus  power-law  (Ramberg-Osgood)  with  an  elliptic 
transition  region.  The  linear  segment  is  necessary  in  the  modified  boimdary 
layer  to  prevent  boundary  plasticity.  The  elliptic  trarisition  region  provides 
continuity  through  the  first  derivative.  The  constitutive  model  is  as  follows 
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(2.15) 

where  coordinates  of  the  center  of  the  ellipse  representing  the 

transitional  region  in  the  effective  stress/strain  space  and  (Ri,R2)  are  the 

ellipse's  major  and  minor  axis  as  shown  in  Fig.  2.1.  These  four  degrees  of 
freedom  are  determined  by  forcing  continuity  between  segments  up  to  the 
first  derivative  in  the  effective  stress/strain  space  for  fixed  values  of  (kj.kj). 
The  relation  for  the  transition  region  is  derived  in  a  marmer  similar  to  that 
given  in  the  previous  section  for  the  Ramberg-Osgood  material  except  the 
equivalent  effective  plastic  strain  is  derived  from  the  equation  of  an  ellipse 
given  by 


(2.16) 


from  which  the  total  uniaxial  strain  is  obtained  as 


Linear 

Elliptic 
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The  plastic  strain  is  obtained  by  subtracting  off  the  elastic  portion  of  the  strain 


assuming  the  multiaxial  plastic  strain  behaves  similar  in  nature  to  the  uniaxial 
relation,  the  following  is  obtained  for  the  effective  plastic  strain 


Substitution  of  the  previous  expression  into  the  deviatoric  stress /strain 
relation  derived  in  the  preceding  section  based  on  the  effective  plastic  strain 
results  in  the  deviatoric  stress /strain  relation  for  the  elliptic  transition  region 
as  given  above. 

The  displacement-based  theory  of  finite  elements  as  the  name  implies 
results  in  the  determination  of  the  displacements  from  which  the  strains  are 
readily  available.  Given  the  strains,  the  stresses  must  be  determined  from  the 
above  constitutive  equation.  Since  the  elliptic  and  the  power-law  segments  of 
the  deviatoric  stress/strain  relations  are  nonlinear  with  respect  to  the 
deviatoric  stress  they  must  be  solved  numerically. 
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The  solution  of  the  constitutive  equations  is  cast  into  the  equivalent 
effective  stress/strain  space  by  taking  the  individual  inner  products  as 
follows 


^Tcul  ^Toid 

3  Bq  Bq 


9^  2  Co  Co 


(2.20) 


where  is  the  expression  in  the  square  brackets  in  the  constitutive  equations. 
The  equivalent  effective  deviatoric  strain  is  obtained  from  the  above 
expression  and  is  given  by 


£o 


=  -p^ 
3  c^o 


(2.21) 


where  Og  is  as  previously  defined  and  the  effective  deviatoric  strain  is  given 
by 

B 


The  equivalent  effective  deviatoric  stress/strain  relations  are  given  by 


(2.23) 


The  linear  equation  is  solved  directly  and  the  nonlinear  equations  are  solved 
using  Newton's  method.  The  above  equations  are  rearranged  to  obtain 
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with  the  Jacobian  given  by 


(2.25) 


Once  the  equivalent  stress  has  been  determined,  the  deviatoric  stress  can  be 
calculated  by  inverting  the  deviatoric  form  of  the  constitutive  equations 


(2.26) 


and  finally,  the  total  stress  is  calculated  using 


(7..  1  er”* 
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(2.27) 
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2.1.2.I  Material  Stiffness 


The  material  stiffness  is  calculated  by  differentiating  the  deviatoric 
stresses  with  respect  to  each  of  the  deviatoric  strains.  It  should  be  noted  that 
from  the  above  expression,  the  following  relation  may  be  obtained 


(2.28) 


Substituting  this  expression  for  beta  and  differentiating  the  deviatoric  stress 
with  respect  to  the  deviatoric  strain  results  in 
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Using  the  effective  stress/ strain  relations,  the  following  may  be  obtained 

da,  =  yde'^°^  (2.30) 

where 
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Differentiating  the  effective  deviatoric  strain  results  in 
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(2.32) 
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Substitution  of  the  previous  results  into  the  stiffness  relation  gives 
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The  material  stiffness  is  updated  incrementally  in  the  user  material 
subroutine. 


2.1.22  Strain  Energy  Density 

The  strain  energy  density  is  computed  for  use  in  the  /-integral 
calculations.  The  total  strain  energy  density  is  obtained  by  integrating  the 
increment  of  work  per  unit  volume  produced  during  deformation. 

Strain  Energy  Density-SED= J  G.jdeJP^‘^K  (2.36) 

This  may  be  rewritten  in  deviatoric  form  as 

where  p  is  the  hydrostatic  pressure  as  defined  previously.  The  deviatoric 
stress/ strain  relations  may  be  substituted  into  the  above  expression  to  obtain 


1  o 

n,(9)|J.n^(9)|f»4n,(9)|‘^;»+Q,(?)| 


where  the  integrated  quantities  are  given  by 


n  +  1 
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The  constitutive  model,  material  stiffness  and  the  strain  energy  density  have 
be  coded  in  FORTRAN  user  subroutine  called  UMAT.  The  user  subroutine  is 
included  in  Appendix  A. 

2.2  Asymptotic  Analysis 

The  following  section  contains  the  procedures  for  the  analytical 
solution  for  the  potential,  stress,  strain,  and  displacement  fields  for  a  crack  in 
an  elastic-plastic  material.  The  governing  field  equations  are  derived  and 
then  the  solution  technique  is  outlined. 

2.2.1  Derivation  of  Governing  Field  Equation 

A  general  power  series  expansion  was  chosen  for  the  analytical 
representation  of  the  near  tip  fields  in  the  present  work.  The  coordinate 
system  of  choice  is  of  cylindrical  type  and  is  centered  at  the  crack-tip  as 
shown  in  Fig.  2.2.  In  the  equations  that  follow,  we  nondimensionalize  the 
radial  distance  from  the  aack-tip  by  means  of  the  /-integral  and  material  flow 
properties: 


r  = 


ggogpr 

J 


(2.40) 


where  f  is  the  dimensional  radial  distance  and  r  is  the  normalized  distance. 
Note  that  /  is  the  actual  dimensional  quantity — ^we  do  not  nondimensionalize 
/  in  this  derivation  unless  otherwise  stated.  Also,  the  stress  potential,  stress, 
strain,  displacements,  and  strain  energy  density  are  normalized  as  follows: 
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where  Oq  is  the  flow  stress,  Eq  is  the  yield  strain  (Cq/^)  and  a  is  a  material 
constant  that  appears  in  the  Ramberg-Osgood  stress  strain  constitutive  law, 
Eq.  (2.14).  Note  that  the  dimensional  stresses  and  actual  strains  are  barred; 
unbarred  stresses  and  strains  are  nondimensional  imless  otherwise  specified. 

The  stress  potential  is  expanded  in  powers  of  the  radial  coordinate  as 
follows: 


0(r, e)  =  (2.42) 

m=J 

where  Km  is  the  unknown  amplitude  of  the  mill  term,  (^"^{6)  is  a 
dimensionless  angular  function  to  be  determined,  and  Sm  +2  is  the  unknown 
exponent  of  r .  The  term  s,n  can  be  thought  of  as  a  pseudo-eigenvalue  in  the 
sense  that  the  formulation  of  this  problem  does  not  result  in  the  traditional 
eigenvalue  type.  Here  each  term  of  the  series  is  assumed  to  be  separable  in  r 
and  0,  but  when  combined  in  their  totality,  the  resulting  expression  is  no 
longer  separable  in  the  two  in-plane  coordinates.  At  this  point  it  should  also 
be  mentioned  that  no  formal  proof  exists  establishing  the  completeness  of  this 
set  of  functions.  In  the  case  that  the  functions  are  complete,  the  series 
expansion  represents  the  exact  stress  potential.  Conversely,  if  the  functions 
are  not  complete,  the  series  is  an  approximation  for  the  stress  potential. 

The  stress  components  are  derived  to  satisfy  the  polar  form  of  the 
equilibrium  equations 


29 


SCr.{r,8)  \d<yjr,e)  ^  ^ 

dr  r  do  r 

^  I  ^°re(''-^)  ^  0 

dr,  r  do  r 

for  the  plane  problem  from  the  potential  as 
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Substituting  the  power  series  expansion  of  the  stress  potential  into  Eq.  (2.44) 
results  in  the  following  stresses 

r’"  (2.45) 

m=i  " 

where 

=  (•'-«+2)(»m  +  l)l’m(9)  <?.46) 

<9(e)  = 

The  expressions  obtained  for  the  stress  field  are  then  used  to  determine 
the  Mises  or  effective  stress 


1 

^Ar.e)  =  [|S(,(r,e)  S..(r.e)p  (2.47) 

where  Sjy(r,0)are  the  deviatoric  stresses  and  are  defined  as  follows 
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Sj,.(r,e)  =  C^(r,e)  -  |oa(r,0)5i,..  (2.48) 

Here  repeated  indices  are  understood  to  be  summed  and  Sy  is  the  well  known 
Kronecker  delta  whose  value  is  =  l  for  i  =  j  and  5y  =  0  when  i  ^  j.  After 

substitution  of  the  series  representation  of  the  stress  fields,  the  deviatoric 
stresses  may  be  written  as 


where 


^<-,6)  =  wm  /”■ 

m=l  J 
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(2.49) 


(2.50) 


The  deviatoric  stresses  are  dependent  on  the  out-of-plane  normal  stress 
a^(r,0)  at  this  point  so  the  third  constitutive  equation  is  solved  to  obtain  the 

following  relation  for  the  out-of-plane  normal  stress 
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where 
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(2.51) 


(2.52) 


Here  lambda  is  boimded  by  v<X(r,0)<O.5.  At  this  point  it  is  worth 
examining  the  situations  for  which  the  limiting  values  of  X  are  approached. 
First,  for  the  case  of  an  elastic  material,  a  -»  0  and  X  v.  Second,  when  a 
material  is  incompressible,  v  ->  0.5  and  X  ->  0.5.  Finally,  when  v  «  a”“  as 
in  a  fully  developed  plastic  region,  X  0.5.  The  use  of  the  above  expression 
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for  the  out-of-plane  normal  stress  in  the  remaining  analysis  results  in  the 
ability  to  solve  both  the  plane  stress  as  well  as  the  plane  strain  problem.  In 
both  cases  is  taken  as  X  times  the  sum  of  the  in-plane  normal  stresses.  For 
the  case  of  plane  stress  we  set  X  equal  to  zero  and  for  plane  strain  we  assume 
X  is  0.5.  The  latter  is  an  assumption  which  makes  the  analysis  somewhat 
neater,  but  will  be  checked  after  the  full  field  results  are  obtained. 

Introducing  the  deviatoric  stresses  in  the  effective  stress,  the  following 
series  representation  is  obtained 

of{r,e) 


where 

(0)  =  fs^(e)  s^{e).  (2.54) 

im  ^ 

Note  that  the  above  expression  represents  the  coupling  angular  fxmctions  for 
the  effective  stress  and  have  units  of  stress  squared. 

The  next  step  in  the  analysis  is  to  express  the  strains  in  terms  of  a 
power  series  by  substituting  the  stresses  into  the  constitutive  relation.  First, 
before  this  can  be  accomplished,  the  effective  stress  must  be  raised  to  the 
power  (n-7).  Application  of  the  binomial  expansion 


(2.55) 

results  in  the  following  series  representation  for  the  expanded  effective  stress 
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where  ^re  angular  functions  dependent  on  a  {6).  The  strains  are 

ij 

obtained  by  introducing  the  above  series  expansions  of  the  stress  components 
into  the  constitutive  equation.  The  series  representation  of  the  strains 
becomes 


(2.57) 


where  E  indicates  the  term  results  horn  an  elastic  field  and  P  from  a  plastic 
field.  are  fimctions  of  and  the  hardening  exponent  n  .  The  exponents 

p  p 

’  are  functions  of  the  stress  exponents  and  the  hardening  exponent  in 
the  plastic  case. 

The  displacement  fields  are  readily  available  from  the 
strain/displacement  relations 


'ee 


r  r  do 


(2.58) 


1  du^{r,e)  du^{r,e)  u^(r,e) 


do 


dr 


The  first  relation  can  be  integrated  directly  to  obtain  the  expression  for  the 
radial  displacement.  The  expression  for  the  shearing  strain  can  be  rearranged 
by  miiltiplying  by  two  and  dividing  by  r  such  that  the  terms  containing  the 
angular  displacement  may  be  combined  under  one  common  partial  with 
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respect  to  the  radial  coordinate,  — The  resulting  expression  can  then  be 

dr 

integrated  with  respect  to  r  to  obtain  the  angular  displacement.  The  resulting 
expressions  for  the  displacement  field  are 

Ur{r,e)  =  ^err{^,e)dr 

rlF  1  (2-59) 

Us(r,0)  =  rj-  2e^^(r,e)-i— dr. 


The  above  series  representation  of  the  strain  field  can  be  substituted  into  the 
integrals  to  obtain 


u,(r,e)  =  *’ +  X^,ufW{e)r'' 


(2.60) 


where 
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The  above  expressions  for  the  displacement  field  ignore  rigid  body 
displacements. 

The  final  equation  to  be  satisfied  in  the  field  problem  is  the  polar  form 
of  the  compatibility  equation 
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Substituting  the  series  expansion  of  the  strains  into  the  Eq.  (2.62)  results  in  the 
following  second  order  power  series  differential  equation: 

f  K,nf{e)r^'^  +  i  K,nf{6y‘'^  =  O  (2.63) 

£=1  £=1 

where 
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The  second  order  differential  equation  is  satisfied  asymptotically  (term 
by  term)  by  solving  the  individual  eigenvalue  problems  =  0.  Each  of 

these  second  order  equations  is  transformed  into  a  fourth  order  ordinary 
differential  equation  in  terms  of  (})'"(6): 


+c;"(e)^^+cj'r(e)+c7(0)=o 


where 


C"‘(0)  =  L7(e)  +  Nf(0) 


(2.66) 


and  the  coefficients  are 


N"(e)  =  N"{n,s„,(|)"(e).(ti”(ey,(l)”'(e)f,(|i”'(e)f"}  (nonlinear ) 
L7{e)  =  L”{n,s„,s^.f  (e),(|i«(ej,4)'(ej',(|i'(e)"'}  (linear ). 

(2.67) 

C”(0)  =  C”{n,s^,()^(0),(i^(0)',4t^(0)",())^{0)"'|  is  the  particular  part  of  the 
differential  equation  and  depends  on  the  solution  to  previous  fields  as 
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designated  by  the  indices  i.  The  differential  equation  corresponding  to  the 
first  term  of  the  series  (HRR  field)  is  homogeneous  in  (j)^(0)  and  and  is 
nonlinear.  The  fields  that  follow  are  inhomogeneous  (i.e.  include  coefficients 
that  depend  on  solutions  to  previous  fields)  but  are  linear  (i.e. 
Nj"(0)  =  O  for  m>l). 

The  solution  to  the  above  differential  equation  is  subject  to  the 
following  boimdary  conditions.  Loading  of  the  Crack  in  mode  I  results  in  the 
following  four  conditions.  The  first  two  are  symmetry  conditions  on  the 
crack  plane.  Since  the  stress  state  on  either  side  of  the  crack  plane  must  be 
identical,  the  stress  potential  must  be  an  even  function  about  0  =  Oand 
therefore  odd  derivatives  of  the  potential  are  odd  fimctions  as  given  by 


5'4>(r,+e)  _  ,  5'3>(r -e) 

50'  ^  '  50' 


where  ^  =  0, 1, 2, 3, . . . .  This  results  in  the  conditions  that 
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The  second  set  of  boundary  conditions  result  from  the  stress  free  surface  on 
the  crack  flank 
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Ore{r,7t)  =  0 


^^0(r.;r)_Q 

54^^ 

d9 


(2.70) 


The  above  conditions  on  the  stress  potential  put  the  following  restrictions  on 
the  angular  functions  in  its  series  representation: 
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de 

d^e 


0  aS(0)  =  0 

0  symmetry 


(2.71) 


Q-m 


d^”^{7t) 

de 


=  0  ->  a2‘g(;r)  =  0 
=  0  ^  c^^iK)  =  0. 


(2.72) 


Combining  these  conditions  with  the  above  differential  equations 
reduces  the  problem  to  determining  the  pseudo  eigenvalues  and 
corresponding  angular  functions  (pseudo  eigen  functions)  which  satisfy  the 
above  boimdary  conditions.  The  solution  of  this  problem  leaves  only  the 
coefficients  in  the  series  expansion  of  the  stress  potential  to  be  determined. 


2.2.2  Numerical  Solution  of  Governing  Field  Equation 


The  fourth  order  differential  equation  derived  in  the  previous  section 
(Eq.  2.65)  was  solved  using  a  multiple  shooting  technique.  The  first  step  in 
the  solution  of  the  differential  equation  was  to  write  it  as  a  set  of  four 
equivalent  first  order  equations  coupled  with  a  single  first  order  equation  for 
the  eigenvalue  as  follows 


37 


d6^ 


VW  =  Yp) 

^  =  Y'm 
d9  ~  de  ’ 

M  =  ^  =  y^(e) 
de^  ds  ^ 


rf¥(£l  =  =  y‘(9\ 

de^  </e  ' 


de 


c'(e)Y'(e)+c'(0)Y'(e)+c'{e)Y'{e)' 

+cJ(e)Y'(e)+c5,(e) _ _ 


c'(e) 


(2.73a) 

and 


S,  =  Y'(e) 

^  ffYMf)  ^  0 
de  de 


(2.73b) 


These  equations  were  then  integrated  from  the  left  (0  =  0)  and  right  (0  =  ^) 
boundaries  using  DIVPRK  a  sixth  order  explicit  Runge-Kutta-Verner 
integration  method  in  the  IMSL  MATH/LIBRARY  [40].  The  solution  at  6  +  h 
is  estimated  by 


rU0+k)=rHo)+ 


— k^’^+Ok^’^ 

40  '  ‘ 


875  , 

- k 

2244  ‘ 


IS 


125 


1955  ‘  '  ■  11592 


23 

+— k 

72 
ky  H 


£.4 


— k^-*l 
160  '  J 


(2.74) 


where  the  indices  £  is  the  field  number,  i  is  the  first  order  equation  number. 
The  coefficients  kp  are  defined  in  Appendix  B.  The  value  h  in  these 
constants  is  the  step  size  for  integration.  An  irutial  step  size  of  0.01  was 
selected  and  was  adjusted  between  0.000001  <  h  <  0.2  to  keep  the  global  error 
<0.00001. 
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The  boundary  conditions  imposed  on  the  above  system  are  shown  in  the 
Fig.  2.3.  The  homogeneous  boimdary  conditions  are  a  result  of  the  symmetry 
and  stress  free  boundary  discussed  above.  The  magnitude  of  the  angular 
fields  in  the  stress  potential  were  arbitrary  set  to  one  at  0  =  0  for  the  purpose 
of  integration  but  were  later  scaled  such  that  the  corresponding  angular  field 
associated  with  the  equivalent  stress  would  have  a  maximum  value  of  one.  It 
was  possible  to  do  the  scaling  because  the  coefficients  in  the  series  expansion 
are  yet  to  be  determined  and  will  scale  the  fields  to  their  appropriate 
magrutudes.  The  boundary  conditions  designated  by  gt  were  initially 

unknown.  It  was  possible  to  determine  their  appropriate  values  by  initially 
guessing  values.  The  correct  values  were  then  arrived  at  by  systematically 
adjusting  them  in  successive  iterations  until  all  the  differences  between  the 
corresponding  functions  Yf(0)  integrated  from  0  =  0  and  those  integrated 
from  d  =  n  had  differences  defined  as 

F|(gj)=  Um  [y/(0)]-  lim  [Yf(0)]  (2.75) 

fi  -4  jr  /  2”  B-^n  12^ 


that  were  less  than  0.0000001.  The  adjustment  of  the  five  unknown  boundary 
conditions  were  determined  through  a  first  order  Taylor  series  expansion  of 
the  difference  functions  F,^(9y)  as  follows 


(2.76) 


where  the  indices  n  and  n+1  indicate  the  current  and  next  iterations, 
respectfully.  The  Jacobian  of  F^(g*.)  was  evaluated  nximerically  for  each 

p 

iteration  by  varying  gj^  individually  by  an  increment  of  Ag^  =0.000001  5^ 
(m  =  1,2, 3, 4, 5)  to  determine  +  Agj(,j.  After  obtaining  these  values, 

the  Jacobian  was  estimated  as  follows 


FIG.  2.3 —  Schematic  of  multiple  shooting  technique. 
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(2.77) 


f 

The  increments  of  the  unknown  boundary  conditions  Ag^were  obtained 

using  the  IMSL  subroutine  DLINRG.  This  subroutine  computes  the  solution 
of  a  linear  system  using  the  LU  decomposition  method.  First,  an  upper  and 
lower  triangular  matrix  whose  product  returns  the  coefficient  matrix  are 
determined  using  Grout's  algorithm. 

Ax  =  LUx  =  b.  (2.78) 


The  solution  of  the  system  is  then  obtained  by  forward  substitution  of  the 
lower  triangular  matrix  and  the  right  hand  side  vector  to  obtain  an 
intermediate  vector. 


Ux  =  L"^b  =  c.  (2.79) 

Next,  back  substitution  of  the  upper  triangular  matrix  and  the  intermediate 
vector  is  performed. 


x  =  U-^c  =  U-^L-^b.  (2.80) 

As  stated  before,  the  differential  equations  were  integrated  successively  each 

f 

time  adjusting  the  values  of  gf  until  the  differences  in  the  right  and  left 

hmctions  Yj(6)  evaluated  at  6  =  7tl2  were  less  than  0.0000001. 

The  angular  functions  for  the  potential  were  then  used  to  determine  the 
angular  functions  for  the  stress,  strain  and  displacement  fields  using  the 
relations  given  earlier.  The  fields  were  considered  correct  if  the  /-integral  was 
path  independent.  The  programs  and  input  files  developed  to  solve  the 
differential  equations  are  included  in  Appendix  C. 
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2.3  Finite  Element  Analysis 

Crack  tip  fields  were  obtained  for  the  modified  boundary  layer  and 
finite  geometry  specimens  using  the  commercial  finite  element  program 
ABAQUS  5.2  [41].  Specimen  modeling  and  post  processing  of  the  results  was 
accomplished  using  the  solid  modeler  PATRAN  2.5  [42].  Data  was  translated 
between  the  two  programs  using  PATABA  and  ABAPAT  versions  3.1A. 

2.3.1  Solution  Methods 


The  field  problem  was  solved  using  the  procedures  of  the  method  of 
virtual  work  rate.  The  principle  is  derived  from  equilibrium  as  follows:  First, 
the  partial  differential  equations  of  equilibrium  in  the  deformed  configuration 
are  multiplied  by  a  test  function  (virtual  velocity) 


■5  +  fl-5v  =  0. 


(2.81) 


In  this  section  arrows  are  used  over  tensors  to  indicate  their  order.  Second, 
the  residual  is  integrated  over  the  domain  of  interest 


m- 


+f\dvdV  =  0. 


(2.82) 


Third,  the  product  rule  is  used  to  rewrite  the  first  term 

Fourth,  the  divergence  theorem  is  used  on  the  first  term  to  obtain 


(2.83) 


-  ddw 


dV-j[T-n5v]dS-j  [f-^]6?V  =  0.  (2.84) 
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Finally,  the  gradient  of  the  virtual  velocity  is  related  to  the  virtual  rate  of 
deformation  by 

^  =  =  +  (2.85) 

where 

SD  and  <5W  .  ■  (2.86) 

sym  2  2 

Here  5D  is  the  rate  of  deformation  tensor  and  5W  is  the  rate  of  spin  tensor. 

Substituting  these  results  into  the  above  expression  and  taking  the 
Cauchy  stress  tensor  x  to  be  symmetric  and  thus  when  multiplied  by  a  skew 
matrix  5W  results  in  the  zero  matrix  the  weak  form  of  equilibrium  is 
obtained 

{x  SDdV-  ft-5vdS-  ff-5vc?V  =  0.  (2.87) 

JV  JS  Jv 

The  statement  of  the  principle  is  as  follows; 

The  rate  at  which  the  external  applied  loads  do  work 
must  balance  the  rate  at  which  the  internal  forces  do 
work.. 

The  first  integral  in  the  above  expression  contains  the  conjugate  work  pair  of 
stress  and  strain  rate.  The  integral  can  be  transformed  so  its  evaluation  may 
be  taken  over  the  original  configuration  provided  the  appropriate  conjugate 
work  pair  are  used.  The  conjugate  pairs  of  interest  are  as  follows 

Stress  Strain  Rate  Domain 

Cauchy,  x  Deformation  Rate  Tensor,  D  current 

Kirchhoff,  d  Deformation  Rate  Tensor,  D  original 

2Dd  piola-Kirchhoff,  S  Green -Lagrange  Strain  Rate,  i  original 

Since  the  rate  behavior  of  materials  is  generally  modeled  with  constitutive 
relations  using  the  deformation  rate  terxsor,  and  we  expect  the  elastic  behavior 
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to  be  derivable  from  a  thermodynamic  potential  defined  over  the  original 
configuration,  the  Kirchhoff  stress  (d)  and  the  deformation  rate  tensor  (5) 
are  used  in  ABAQUS.  The  Kirchhoff  stress  is  related  to  the  Cauchy  stress  (  t) 
by 


?  = 


g 

detlS 


(2.88) 


-  dx 

where  F  is  the  deformation  gradient  and  x  are  the  deformed  coordinates 


while  X  are  the  tmdeformed  coordinates.  Although  the  Kirchhoff  stress  is 
used  in  the  constitutive  model,  the  final  stress  output  is  the  Cauchy  stress. 

The  finite  element  form  of  the  principle  of  virtual  work  rate  is  obtained 
by  discretizing  the  original  domain  into  a  finite  number  of  subdomains  called 
elements.  The  displacement  formulation  interpolates  the  displacement  field 
within  each  element  from  locations  (nodes)  where  the  displacements  are 
sought.  The  principle  of  the  virtual  work  rate  is  applied  over  each  element  in 
the  following  form 


fBMddV„-fNj,IdS-fN;[,-f<A^  =  0  (2.89) 

vVp  JS  •'V 


where  Nj^  are  the  interpolation  (shape  functions)  and  is  the  displacement 
to  strain  transformation  matrix  containing  spatial  derivatives  of  the 
interpolation  functions.  The  contributions  of  all  elements  are  then  summed. 
The  above  equation  is  in  general  nonlinear  in  nature  and  can  be  written  as 

F'^(u^)  =  0  (2.90) 

« 

where  are  the  unknown  displacements.  The  nonlinearities  result  from 
four  sources;  material  nonlinearity,  large  displacement  gradients,  large 
displacements  and  rotations,  and  boundary  conditions.  The  nonlinear  system 
is  solved  incrementally.  First,  the  equation  is  expanded  in  a  Taylor’s  series 
about  the  approximate  solution  uf  as 
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The  second  term  contains  the  Jacobian  (stiffness)  matrix  for  the  nonlinear 
system  and  is  given  by 


Kmn  =  IBm^ Bn-^Vo  +  -  JnJ,  Q^dS- JnJ, 

(2.92) 

where  H  is  the  current  material  tangent  stiffness  tensor.  The  integrals 
containing  and  result  from  what  are  termed  follower  forces  and  are 
generally  called  the  load  stiffness  matrix  and  are  defined  as 


os  -  JL+_Li;^ 

qV  -  ^ 

J<?Un 


(2.93) 


where  is  the  ratio  of  surface  areas  between  the  current  and  reference 
configurations  and  J  is  the  ratio  of  the  volume  between  the  current  and 
reference  configurations. 

Successive  corrections  are  made  to  the  nodal  displacements  by  solving 
the  linear  system 


After  each  iteration,  the  updated  displacement  vector  is  used  to  determine  the 
residual  in  the  nonlinear  system.  When  the  residual  is  below  some 
predetermined  convergence  tolerance,  the  system  has  been  solved  for  that 
increment. 

After  convergence  is  obtained,  the  strain  increment  is  calculated  from 
the  rate  of  deformation  using  the  following  relation 
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f-e  ^t+At 

Ae  =  e-e^  =  \de=  Ddt 
^0 

and  is  related  to  the  velocity  gradient  through 

where  L  is  the  velocity  gradient  defined  over  the  increment  of  time  At  as 

L(t+T)  =  AR(t  +  At)[AR(t  +  T)'^  L  AR(t  +  T)]AR(t+At)^.  (2.97) 

Using  the  polar  decomposition  of  the  deformation  gradient  increment,  the 
part  in  brackets  in  the  previous  expression  once  substituted  in  the  integral  can 
be  rewritten  in  terms  of  the  increment  of  the  right  stretch  tensor 

ARf  (L + L^jAR  =  AUAU-'  +  AU-'AU  .  (2.98) 

Now  assuming  that  the  incremental  stretch  has  the  same  principle  directions 
throughout  the  increment  we  can  write  the  increment  of  stretch  as 

AC/ = fl + {AX,  -  .  N,  +  (l + {AX„  -  .  N„ 

+ ^1 + {AXjjj  -  •  ^ur 

Substituting  into  the  integral  and  integrating  results  in  the  logarithmic  strain 
increment  given  by 

Ac  —  In^AXfjjiij '  n^  *!■  In^AXjj^njj  •  Hjj  4"  In^AX'jjj  Jnjjj  ■  Hjjj  (2.100) 

where  X,-  are  the  eigenvalues  of  the  left  stretch  tensor  F  •  and  N,-.  n^-  are 
the  eigen  vectors  at  the  beginning  and  end  of  the  increment,  respectfully  (i.e. 
different  by  AR).  The  strain  increment  is  sent  along  with  all  necessary 
variables  to  the  constitutive  routine  to  obtain  the  increment  of  Kirchhoff 


(2.95) 


(2.96) 
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stress.  This  stress  is  ther\  transformed  to  the  Cauchy  stress  in  the  deformed 
configuration. 

All  specimens  were  modeled  using  continuum  plane  strain  eight  node 
isoparametric  elements  (CPE8).  These  elements  are  displacement  based  and 
use  the  same  order  shape  functions  for  both  variable  interpolation  within  the 
element  and  coordinate  mapping  between  the  local  (master)  and  global 
spaces.  The  mapping  of  coordinates  and  interpolation  of  displacements  are 
given  as  follows 


^  (2.101) 

where  j  is  an  indices  over  the  number  of  nodes  in  the  element,  xf  are  the 
global  coordinates  of  the  nodes,  u/  are  the  nodal  displacements,  N^.  are  the 
shape  functions  defined  over  the  master  element  given  in  Appendix  D,  are 
the  local  coordinates,  and  x,-  and  are  the  global  coordinates  and 
displacements,  respectfully.  A  3  x  3  Gauss  quadrature  with  nine  integration 
points  was  used  for  integrating  over  the  master  element.  Full  integration  was 
deemed  necessary  to  capture  the  complete  effect  of  plasticity  near  the  crack- 
tip.  Reduced  integration  which  is  used  frequently  to  reduce  the  effect  of  over 
estimating  the  stiffness  matrix  in  the  displacement  formulation  was  not  used. 
The  calculations  where  conducted  with  the  following  expression 


Q  m  =  ;Ln  =  i  J 


(2.102) 


where  are  the  coordinates  of  the  nine  quadrature  points  and  are  the 
nine  quadrature  weights  their  values  are  given  in  Appendix  D.  A  schematic 
of  the  CPE8  element  is  given  in  Fig.  2.4  with  both  its  node  and  integration 
point  numbering. 
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2.3.2  Modified  Boundary  Layer  Models 


A  modified  boundary  layer  as  shown  in  Fig.  2.5  with  a  displacement 
field  corresponding  to  the  first  three  Williams  terms  {Kj  T  and  S)  was 
constructed  using 


Uj,{f,e)=-^Y  cos^|j{3-4v-cos(0)} 


1  1  T 
f2+. 


(l-v^) 


cos{d)r 


+  5 


K,{l+v) 

Uy(r,e)=-^ - 


(1  +  v) 


cos(0)n --4v Jcos|^-j j) 


(2.103a) 


E 


+  5 


in^|j{3-4v-cos(a)}  f2  ^sin(0)f 

sin(  -  4 '' j  ^ 

+cos(e)([i-4v)sin[|]-isin[^|j^ 


(1+v) 


(59' 

y  2  . 


(2.103b) 


The  first  two  terms  (X/  and  T)of  the  Williams  solution  was  chosen  as  our  two- 
parameter  field,  Larsson  and  Carlsson  [43].  All  higher  order  coefficients 
where  set  to  zero.  The  evolution  of  the  plastic  field  in  this  model  by 
definition  depends  only  on  two  parameters  Kj  and  T  ,  excluding  material 
properties.  For  this  reason,  the  second  and  third  independent  coefficients  in 
the  nonlinear  asymptotic  analysis  must  only  depend  on  the  same  two 
parameters  Kj  and  T  when  the  analytical  representation  for  the  stress  field  is 
matched  with  the  full  field  finite  element  results.  A  separate  model  was 
constructed  to  study  the  influence  of  the  third  Williams  term  on  the  stress 
field  near  the  crack-tip.  Displacements  corresponding  to  Kj  and  S  where 
applied  to  the  boimdary  in  this  model. 

The  region  was  modeled  as  a  25.4  cm  semi-circle,  2.54  cm  thick.  There 
are  54  rings  of  elements  which  geometrically  degenerate  towards  the  crack- 
tip.  The  ratio  of  the  radial  dimension  of  the  outermost  element  ring  to  the 
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inner  most  ring  is  3.21251  cm  /  0.001145  cm  =  2805.69.  There  are  16  rays  of 
elements  centered  at  the  crack-tip  with  an  angular  measure  of  0.19634 
radians.  The  aack-tip  elements  were  degenerated  to  triangles.  The  mid-side 
nodes  where  left  in  place  and  the  crack-tip  nodes  were  not  tied,  resulting  in  a 
1  /  r  strain  singularity.  The  analysis  was  based  on  small  strain  theory. 

The  bovmdary  displacements  in  the  first  model  were  applied  by  first 
applying  displacements  associated  with  a  X;  value  of  54.96MPaVm  and  then 
displacements  associated  with  T/ob  varying  from  -1  to  +1  in  0.1  increments.  It 
was  confirmed  by  using  27.48MPa'\/m  as  an  alternate  value  of  X/  that  the 
stress  distribution  depends  only  on  the  value  of  T/ob  when  the  radial  distance 
is  normalized  by  //ob.  Small  scale  )delding  fSSYj  was  maintained  by  insuring 
that  the  value  of  the  numerically  measured  /  was  at  least  99%  of  J  determined 
from  X/ applied  at  the  boundary.  This  was  accomplished  in  two  ways.  First, 
the  plastic  zone  rp  was  restricted  to  ten  percent  of  the  radius  of  the  model  and 
second  a  constitutive  model  was  used  which  behaved  in  a  linear  elastic 
manner  up  to  =  0.950^  and  then  became  nonlinear  so  that  at  large  strains  the 
stress/strain  curve  corresponded  to  the  Ramberg-Osgood  curve. 

The  boundary  displacements  in  the  second  model  were  applied  in  a 
marmer  similar  to  that  of  the  first  model.  The  applied  X/  value  was  reduce  to 
27.48MPaVm  to  maintain  small  scale  yielding.  T/ob  was  set  to  zero  and  S/ob 
was  varied  from  -1  to  +1  in  0.1  increments,  /-integrals  were  calculated  by  the 
domain  integral  method  in  both  models  and  showed  path  independence  for 
contours  taken  outside  the  finite  strain  zone.  A  sample  ABAQUS  input  file  is 
included  in  Appendix  E. 

2.3.3  Finite  Geometry  Specimen  Models 

The  finite  geometry  specimens  are  the  single  edge  notch  bend  (SENB), 
single  edge  notch  tension  (SENT),  center  cracked  tension  (CCT),  and  double 
edge  notch  tension  (DENT).  Figs.  2.6  and  2.7  show  the  mesh  and  boundary 
conditions  for  each  of  these  geometries.  Since  the  three  specimens  could  be 
modeled  using  a  similar  geometry  by  only  changing  the  boundary  conditions, 
the  following  dimensions  apply  to  all  of  the  specimens;  the  depth  of  the 
specimens  (W)  was  5.08  cm,  a  crack  length  to  width  ratio  (a/W)  of  0.1,  0.5, 
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and  0.9  were  used,  the  width  to  thickness  ratio  (W /B)  was  2.0,  and  the  half 
length  (L)  of  the  specimen  was  10.16  cm.  The  focused  mesh  used  in  the 
modified  boimdary  layer  analysis  was  scaled  and  inserted  for  the  crack-tip 
region. 

Loads  were  applied  to  the  specimens  such  that  the  plastic  J  values,  as 
determined  by  the  EPRI  method  [44],  were  given  as  Jpi  =  aoo/A  The 
coefficient  ^  was  chosen  to  be  20  for  all  geometries  except  the  CCT  which  was 
loaded  to  a  P  value  of  50.  The  resulting  loads  were  calculated  with  the 
following  expression 

1 

n+l 

(2.104) 


P=P, 


Pae^hh^l 


where  P©  is  the  specimens  limit  load,  n,  £o,  a,  are  material  constants  from  the 
Ramberg-Osgood  constitutive  relation,  b  is  the  uncracked  ligament  length, 
and  hi  is  the  EPRI  coefficient  for  a  given  geometry.  The  load  was  applied 
incrementally  in  twenty  equal  steps.  In  all  but  the  SENB  specimen  this  load 
was  distributed  evenly  over  the  surface  on  which  it  was  applied.  Stresses 
were  determined  for  a  cylindrical  coordinate  system  centered  at  the  crack-tip. 
/-integrals  were  calculated  by  the  domain  integral  method  and  showed  path 
independence  for  contours  taken  outside  the  finite  strain  zone.  A  sample 
ABAQUS  input  file  is  included  in  Appendix  E. 


2.4  J-Integral  Analysis 


The  /-integral  is  a  well  known  parameter  used  to  characterize 
nonlinear  material  behavior  ahead  of  a  crack,  it  was  first  discovered  as  a 
general  conservation  integral  by  Eshelby  [45]  and  was  then  applied 
independently  to  the  crack  problem  by  Rice  [9].  The  value  of  the  /-integral 
gives  a  measure  of  the  energy  flux  focused  on  the  crack-tip.  It  is  derived  here 
from  an  integral  over  a  two  dimension  area  (shown  in  Fig.  2.8)  containing  no 
singular  behavior  given  by 
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f(WA  =  0. 
Ja 


(2.105) 


This  expression  holds  for  any  area  of  interest.  Substituting  for  the  integrand  a 
seemingly  arbitrary  term  results  in 


f  \Ma  iil— 


(2.106) 


Passing  the  differential  operator  through  the  brackets  on  the  first  term  results 


(2.107) 


Noting  the  first  term  contains  the  equations  of  equilibrium  and  using  the 
symmetry  of  the  stress  tensor  along  with  changing  the  order  of  integration  on 
the  displacements  in  the  second  term  the  expression  can  be  rewritten  as 


f  L  5  \i/^i ,  ^A1  d  f. 


=  0.  (2.108) 


de.. 


The  first  term  may  be  rewritten  as  ^ij~^  which  for  materials  whose 

dW 

stress /strain  relation  can  be  derived  from  a  strain  energy  potential,  =  a^y 
the  entire  expression  can  be  written  as 


r  r dW  d  f-  1 

6..^  dA  =  i 


(2.109) 


Using  the  divergence  theorem  to  convert  the  area  integral  to  a  line  integral 
results  in 
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duj 


k/5  =  0 


and  finally  eliminating  the  direction  cosines  gives 


-t^dS 

idXj 


=  0. 


(2.110) 


(2.111) 


Taking  the  contour  S  to  be  that  shown  in  Fig.  2.8  the  integral  can  be  split  into 
four  separate  integrals  along  Fj-  where  S=r^+r2  +  r_j  +  r^.  Noting  that  the 

paths  along  the  crack  faces  do  not  contribute  because  dy  and  t/  are  both  zero 
results  in 


"I 

Therefore  when  the  contour  is  taken  as  an  open  contour  intersecting  the  crack 
surfaces,  the  expression  for  the  /-integral  is  path  independent  and  can  be 
expressed  as 


Wdx,-t.^dr, 

2  tdxj  ^ 


=  0. 


(2.112) 


/  = 


r 


.  du.  ^ 
-T.^dS 
^  dx. 


(2.113) 


The  /-integral  was  used  in  the  analytic  as  well  as  the  numerical  portions  of  the 
present  work. 


2.4.3  Analytic  J-integral 

The  /-integral  was  used  in  the  analytic  analysis  to  nondimensionalize 
the  radial  coordinate  as  well  as  to  check  the  angular  stress  and  displacement 
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fields.  The  power-law  representations  for  the  stress,  strain  and  displacement 
field  were  substituted  into  the  equation  for  /  given  above  to  obtain  a  general 
expression.  The  result  is  a  series  expansion  of  /  in  which  all  the  terms  but  the 
first  have  exponents  on  the  radial  coordinate  other  then  zero.  For  /  to  be  path 
independent  the  coefficients  of  these  terms  must  be  zero.  Each  of  the 
coefficients  contain  an  integral  with  limits  of  -ii<Q<n  whose  integrand 
contains  the  angular  functions  of  stress,  strain,  and  displacement.  The 
angular  fields  were  considered  correct  if  the  integrals  were  zero.  The 
following  is  the  method  by  which  /  was  determined  for  the  analytic  work. 

The  expression  for  the  /-integral  can  be  expressed  in  polar  form  by  first 
substituting  for  the  traction  vector  and  then  expanding  to  obtain 


[(  3u 

3u,] 

f 

du. 

3u  ] 

^  1 

J  - 

i.a  _ 2 

n  4- 

/T 

1 

4-  ^  ^ 

"2 

^ds 

dx 

^12 

dx. 

IV  1 

i) 

K 

1 

1 ) 

(2.114) 


where  the  term  with  has  been  omitted  because  =  0 .  The  field  variables 
are  transformed  following  the  rules  for  tensor  transformations  as  follows 


’cos0  -sin0' 

«2(r,e) 

sin0  COS0 

u^[f,e\ 

(2.115) 


d„(r,0) 

cos^  6 

sin^  B 

~2sin0cos0 

a-j2(r,e) 

,  — 

sin^  B 

cos^  B 

2  sin  0  cos  0 

>.  (2.116) 

sin  0  cos  0 

~sin0cos0 

cos^  0-sin^  0 

The  gradient  in  the  displacement  field  may  be  expressed  in  polar  form  vising 
the  chain  rule 


where 


dx.  do  dx  dr  dx 

^  ^  ^  (2.117) 

du^(r,e)  _  X^(r,0]  30  du^lKd)  dr 

dxj  ~  de  dXj  3r  dxj 


0  =  tan  ^ 

f  \ 

X 

2 

A 

X, 

/i/Tv 

(2.118) 

L  ij 

do  _ 

— sin0 

dr 

A 

r 

* 

=  COS0 

(2.119) 

1 

II 

II 

COS0  ds. 

(2.120) 

The  polar  form  of  the  /-integral  for  a  circular  contour  is  obtained  by 
substituting  the  above  expressions  into  the  Cartesian  form  of  the  integral. 
After  simplification,  we  obtain 


. .  =  f  W(r,e)cosd-cosd  a 

aa  e  J  '  '  '  dr 

0  0  ^  L 


-cose 

-2|»a„(f.e)j«,(f,e)-i^ 

de 


(2.121) 


fde 


where  the  field  variables  are  the  dimensional  values.  The  strain  energy 
density  was  obtained  by  integrating  for  the  area  under  each  of  the 
stress /strain  relations  and  summing  the  results  as  follows 
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W(r,0)  =  a  (f,e)e..(f,0)-  f 

•'O 


o..(r.  6) 

y 


i..(r.e)dc.(r,e).  (2.122) 


The  strain  energy  density  for  the  Ramberg-Osgood  constitutive  model  can  be 
expressed  as  follows 


(2.123) 


Substituting  the  series  expansion  of  the  stress,  strain,  and  displacement  fields 
into  the  strain  energy  density  and  then  into  the  polar  expression  for  the  f- 
integral  one  obtains  the  following  relation 


T  As  As  As  As  As  As 

—^  =  Ir  1+lP  2^Ir  ^+/r  ^+Lr  ^+/,f  (2.124) 

aOnSr  1  2  3  4  5  6 


’O'^O 


where  As^  ^s^-Sj  and  1/  are  constants  shown  in  Appendix  F.  The  magrutude 
of  the  exponents  increase  with  the  first  exponent  being  zero.  Hence,  as  stated 
above,  for  /  to  be  path  independent  (i.e.  for  any  value  of  f ),  the  constants  li 
must  be  zero  for  i  >1.  The  angular  fields  were  checked  by  evaluating  each  of 
the  constants  Ij . 


2.4.2  Numerical  J-lntegral 


The  /-integral  was  evaluated  for  the  modified  boundary  layer  and 
finite  body  specimens  using  the  domain  integral  techivique.  This  method  is 
implemented  within  ABAQUS.  The  domain  integral  form  of  the  /-integral  is 
obtained  by  multiplying  the  integrand  by  a  function  q  as  follows 


(2.125) 
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The  magnitude  of  q(S)  is  dependent  on  the  position  along  the  contour.  It  is 
equal  to  one  on  r_j  and  zero  on  the  rest  of  the  contour.  The  fimction  q  must 

transition  between  the  value  of  one  and  zero  smoothly  over  the  enclosed 
domain.  In  the  general  case,  r_j  must  be  taken  vanishingly  small.  This  is  not 

conducive  to  numerical  algorithms.  For  this  reason,  the  above  expression  is 
transformed  into  an  integral  over  the  enclosed  domain  in  which  q  is  defined. 

The  above  expression  is  converted  by  applying  the  divergence 
theorem.  The  resulting  equation  for  the  /-integral  is 


The  next  step  is  to  use  the  product  rule  on  the  integrand  which  results  in  the 
following 


du. 
a  — ^ 
V  dx. 


U(S)  , 

dw  _ 
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aix  ] 

3x. 

/  1 

’  dx.' 
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) 

LM.  (2.127) 


The  second  term  in  the  integral  is  equal  to  zero  as  shown  previously.  The 
final  expression  for  the  /-integral  is  given  by 

A 

The  integral  is  now  evaluated  numerically  using  Gauss  quadrature.  Each  of 
the  terms  in  the  integrand  is  defined  in  the  finite  element  formulation  by 
interpolation  from  nodal  points  using  element  shape  fimctions.  The  function 
q  is  given  the  value  one  at  nodes  along  F j  and  zero  at  nodes  along  the  rest  of 

the  contour.  The  magnitude  of  q  at  interior  nodes  is  set  at  one. 

Numerical  values  of  the  /-integral  were  obtained  for  twenty  contours 
around  the  crack-tip.  The  first  several  contours  were  not  path  independent 


W5 
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(2.128) 
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which  results  from  high  distortion  of  elements  near  the  tip  of  the  crack  in 
what  is  known  as  the  finite  strain  zone.  After  the  first  six  to  ten  contours,  the 
value  of  /  was  path  independent. 

2.5  Amlytic/Numerical  Field  Matching 


The  coefficients  in  the  general  power  series  representation  of  the  stress 
field  are  determined  by  matching  the  analytic  solution  with  the  full  field  finite 
element  solutions.  The  stress  from  the  analytical  solution  was  matched  on  a 
point  by  point  bases  with  the  modified  boundary  layer  and  finite  body 
results.  The  following  section  contains  an  outline  of  the  procedures  used  to 
determine  the  coefficients. 

The  series  expansion  of  the  stress  field  in  normalized  coordinates  is 
given  by 


The  third,  fifth  and  sixth  coefficients  for  the  values  of  the  hardening 
coefficients  used  in  the  present  study  are  not  independent  but  depend  on 
previous  coefficients  as  will  be  discussed  in  the  results  section.  The  relation 
between  these  coefficients  and  the  previous  are  given  by 

K  K 

^6  ~  iC^'  (2.130) 

Substitution  of  the  above  into  the  equation  for  the  stress  gives 


(2.131) 
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In  the  first  six  terms,  only  three  independent  coefficients  remain.  The  first  is 
determined  from  the  /-integral.  The  remaining  two  must  be  computed  from 
matching  with  the  finite  element  results.  The  expression  for  the  stress  may  be 
rewritten  as 

(r,  e) = 4 + + 4^:1 +^K,+ 4^1 + 4^r,Ar, + ■  ■  •  {2.132) 

where  depends  on  the  first  coefficient,  the  angular  stress  fields,  and  the 

radial  coordinate.  The  solution  of  the  two  remaining  coefficients  requires  an 
independent  equation  for  each.  In  this  work,  the  radial  and  tangential 
stresses  from  the  finite  element  results  were  used.  The  following  equations 
were  obtained 


—FEA(^  o\~.  jl  >  t3  j^2  >j4j^  i  rS  j^3  ^  76  v  v 

+  ^^4  +  ^^2  + 


(2.133) 


The  fourth  coefficient  was  eliminated  from  these  equations  resulting  in  the 
equation  for  K2  as  follows 

+  G,Kl  -I-  G,Kl  +  G,K^  +  G,  =  0  (2.134) 

where 


^2  -  ^bB^  ^BB^  ~  ^BB^r  ~  ^BB^r 
~  ^BB^r  ^BB^  “  ^6B^  ~  ^B6^ 

G.  =  +4,Z^,-4,Zi-4,i« -4<7™(r.e)+Zt(T™(r,e) 

G,  =  4sZ4-4,I^.-4X“('-.®)+i^<^e“(Ge)- 

(2.135) 


K2  was  determined  from  Eq.  (2.134)  using  Newton's  method.  The  FORTRAN 
listing  is  given  in  Appendix  G. 
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CHAPTER  III 
RESULTS 


The  analytic,  numerical  and  comparison  results  are  given  in  the 
following  sections.  First,  the  method  for  determining  the  exponents  in  the 
general  power  series  expansion  of  the  stress  potential  is  reviewed  and  then  the 
exponents  are  presented.  The  angular  fields  for  stress,  strain,  and 
displacement  corresponding  to  the  first  six  exponents  are  introduced. 

Second,  the  full  field  finite  element  results  are  presented.  The  stresses 
obtained  from  the  finite  element  method  for  the  geometries  of  this  study  are 
given  along  the  crack  plane  and  as  a  fimction  of  theta  at  a  normalized  radius 
of  approximately  two.  The  six  term  analytical  representation  of  the  stress  field 
is  also  presented  with  the  numerical  results.  Comparisons  are  made  between 
the  stress  fields  of  the  MBL  for  two  different  applied  K  valves.  The  influence 
of  second  and  third  Williams  term  on  the  stresses  within  the  plastic  region  are 
reviewed.  Next,  the  assumption  made  in  the  third  constitutive  equation  of  the 
analytical  work,  namely  that  X(r,  6)  =  0.5  for  plane  strain  is  addressed. 

Finally,  the  evaluations  of  the  second  and  third  independent 
coefficients  in  the  general  power  series  representation  of  the  stress  field  are 
compared  between  the  MBL  and  the  finite  body  specimens.  These  results  are 
used  to  determine  when  the  stress  field  within  the  plastic  region  of  the  finite 
bodies  may  no  longer  be  characterized  by  a  two-parameter  field. 

3.2  Exponents 

The  angular  fields  and  their  corresponding  exponents  are  obtained  by 
satisfying  the  series  representation  of  the  compatibility  equation  (Eq.  2.63) 
asymptotically  term  by  term.  After  using  the  binomial  exparision  of  the 
effective  stress  in  the  plastic  part  of  the  strain  and  substituting  into  the 
compatibility  equation,  the  exponents  of  the  radial  coordinate  become  integer 
combination  of  preceding  exponents  (i.e.,  =  nsj  +aj.ASj^where  As^ 

and  aj.  are  positive  integers).  The  corresponding  ordinary  differential 
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equation  is  given  by  7r^(0),  Eq.  (2.63).  the  angular  plastic  stress/strain 
relation  is  given  by 


(0)+(«-i) 


(3.1) 


for  terms  vy^here  =  nsj  +  A&-.  The  fourth  order  ODE  (Eq.  2.65)  corresponding 
to  these  terms  is  homogeneous  in  nature.  The  corresponding  exponents  are 
independent  of  any  preceding  fields.  The  terms  associated  with  exponents  of 
the  form  t^.  =  nSj  +  where  aj,  is  greater  than  one  for  single  ASj^  or  k 

takes  on  two  or  more  values  are  not  satisfied  directly  by  the  solution  of  the 
homogeneous  equations.  Hence,  these  terms  are  satisfied  by  grouping 
homogeneous  terms  with  what  are  considered  the  particular  terms.  This 
reordering  of  terms  results  in  the  exponents  being  independent  if  the  value  is 
obtained  from  the  homogeneous  equation  or  dependent  if  the  exponent  on  a 
homogeneous  equation  is  forced  to  be  some  combination  of  previous 
exponents  to  satisfy  a  particular  solution. 

The  exponents  (si)  in  the  general  power  series  expansion  of  the  stress 
potential  are  used  to  satisfy  the  boimdary  conditions  at  0  =  tt,  Eq.  (2.72b).  At 
6  =  0  the  first  and  third  derivative  of  the  angular  potential  fimction  must  be 
identically  zero  to  satisfy  the  symmetry  conditions.  This  leaves  two 
unspecified  boimdary  conditions,  the  value  of  the  angular  potential  field  and 
its  second  derivative,  and  the  exponent.  Either  of  the  two  boundary 
conditions  can  be  fixed  at  0  =  0  to  make  the  problem  well  posed.  Here  we 
chose  to  fix  the  magnitude  of  the  angular  potential  field  to  a  magnitude  of  one. 
Subsequent  to  solving  for  the  field  the  angular  functions  are  rescaled  such  that 
the  value  of  the  angular  equivalent  stress  function  does  not  exceed  unity.  The 
second  derivative  of  the  angular  potential  function  and  the  exponent  are  then 
selected  to  satisfy  the  two  boundary  conditions  at  6  =  n. 

The  first  independent  exponent  is  determined  by  solving  for  the 
leading  nonlinear  angular  field.  This  was  first  solved  by  Hutchinson  [7],  Rice 
and  Rosengren  [8].  The  angular  functions  associated  with  the  first  term  have 
been  coined  the  HRR  field.  The  solution  to  the  first  term  was  obtained  using  a 
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multiple  shooting  technique  as  outlined  in  the  procedures.  The  numerical 
value  obtained  for  the  leading  exponent  was  compared  to  that  derived  by  Rice 
who  used  the  /-integral,  Sj=-l/n  +  l.  The  numerical  values  agreed  with 

Rice’s  expression  up  to  8  sigruficant  figures. 

Once  the  leading  term  was  obtained,  the  higher  order  independent 
exponents  were  determined.  The  differential  equations  governing  the  higher 
order  independent  exponents  are  linear  and  similar  in  nature.  The  fourth 
order  governing  equations  are  obtained  by  substituting  the  plastic  strains  in 
Eq.  (3.1)  into  Eq.  (2.64).  Since  all  higher  order  plastic  strairrs  depend 
exclusively  on  the  present  field  and  the  leading  (HRR)  field  in  the  same 
manner,  the  higher  order  independent  exponents  can  be  determined  from  the 
same  equation. 

The  higher  order  independent  exponents  were  determined  in  the 
following  manner.  Since  the  governing  equation  is  linear,  two  elementary 
solutions  are  sought  which  satisfy  the  following  botmdary  conditions  at  0  =  0 


(1)  <ff{0)=l, 

(2)  (|,”(0)=0, 


d^(|."(0)_ 

’ 


(3.2) 


respectfully.  The  elementary  solutions  for  the  above  conditions  were  obtained 
by  integrating  the  governing  equation  with  the  sixth  order  explicit  Runge- 
Kutta-Verner  integration  method  outlined  in  the  procedure  section  for  a  fixed 
valve  of  the  exponent  The  general  solution  is  then  given  by  combining  the 
elementary  solutions  as  follows 

r(e)=q<l>r(e)+q'l'r(®)  P-s) 

were  Ci  and  C2  are  real  constants.  Boundary  conditions  at  6  =  n  must  also  be 
satisfied.  Hence, 


r(7c)  =  0  =  q({)f(7t)+C2())J(7:) 

de  de  ^2  • 


(3.4) 
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FIG.  3.1 —  Determination  of  independent  eigenvalues  for  plane  strain 
with  a  hardening  coefficient  of  n=5. 
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FIG.  3.2 —  Determination  of  independent  eigenvalues  for  plane  strain 
with  a  hardening  coefficient  of  n=10. 
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These  relations  can  be  rearranged  in  matrix  form  to  obtain 


(j)f(7C) 

.  de 


d9  J 


[01 

>  * 

0 

J 

(3.5) 


A  non-trivial  solution  for  the  vector  {C}  exists  if  and  only  if  the  determinant 
of  the  matrix  of  the  coefficients  vanishes,  i.e. 


de 


(j)”(7C)  =  0. 


(3.6) 


This  is  satisfied  for  discrete  values  of  the  exponent  which  correspond  to  the 
independent  exponents  3.1-3.3  show  Eq.  (3.6)  plotted  as  a  function  of 

the  exponent  value.  Figs.  3.1  and  3.2  are  the  results  for  the  case  of  plane  strain 
with  hardening  exponents  of  5  and  10,  respectfully.  Fig.  3.3  shows  the  results 
for  the  case  of  plane  stress  with  a  hardening  exponent  of  5.  The  locations  were 
the  determinant  passes  through  zero  represent  the  independent  exponents.  In 
the  figures,  the  independent  exponents  are  numbered  by  the  order  in  which 
they  appear  in  the  series.  The  first  cross  over  point  corresponds  to  the  HRR 
exponent  Sj  =  -1  /  n  + 1.  The  first  four  independent  exponents  were  obtained 

for  the  case  of  plane  strain  while  only  the  first  exponent  could  be  determined 
for  the  plane  stress  case.  Xia,  et  al.  [17]  and  Yang,  et  al.  [18,19]  have  reported 
similar  findings  for  exponents  up  to  the  fourth  field. 

The  dependent  exponents  were  determined  by  grouping  terms  in  the 
series  representation  of  the  compatibility  equation.  Some  of  the  homogeneous 
terms  were  grouped  with  those  that  do  not  have  coefficients  that  vanish 
(coefficients  and  exponents  depend  on  solutions  of  earlier  fields  and  hence 
provide  a  particiilar  solution  to  a  homogeneous  equation)  with  independent 
terms.  The  asymptotic  analysis  is  carried  out  by  ordering  the  exponents  on 
the  radial  coordinate  in  increasing  order.  The  ordinary  differential  equation 
that  results  from  the  grouping  is  linear  in  nature  and  is  associated  with  a 
particular  solution.  The  solution  of  these  fields  was  obtained  by  combining 
the  two  elementary  solutions  for  the  homogeneous  equation  with  the  solution 
for  the  particular  part,  i.e. 
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r(0)=q(i)r(^)+cA"M+^pM-  (3-7) 


The  coefficients  were  determined  by  satisfying  the  boundary  conditions  at 
6  =  TT  by  solving  the  system 
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(3.8) 


The  exponents,  their  type  and  ordering  are  listed  in  Tables  3.1-3.3.  Plane 
strain  results  for  hardening  exponents  of  5  and  10  are  given  in  Tables  3.1  and 
3.2,  respectively.  Table  3.3  contains  the  results  for  plane  stress  with  a 
hardening  exponent  of  5. 

3.2  Angular  Fields 

The  corresponding  angular  stress,  strain,  displacement,  and  /-integral 
fields  derived  from  the  first  six  angular  potential  functions  for  plane  strain 
with  hardening  coefficients  of  10  and  5  are  given  in  Figs.  3.4-3.9  and  Figs.  3.10- 
3.15,  respectively.  Each  of  the  angular  fields  were  evaluated  for  errors  by 
checking  for  path  independence  of  the  /-integral.  The  angular  integration 
constants  Ij  as  given  in  Appendix  F  are  shown  in  Figs.  3.9  and  3.15.  The 
constants  in  all  but  the  first  term  integrate  to  zero  which  preserved  the  path 
independence  of  the  /-integral.  The  first  four  angular  stress  fields  agree  with 
those  reported  by  Yang  et  al.  [18,19].  All  but  the  third  field  agree  with  those 
reported  by  Wang  et  al.  [17]  who  recently  have  indicated  there  was  an  error  in 
their  third  angular  stress  fields. 

3.3  Finite  Element  and  Analytic  Comparisons 

The  HRR  or  first  term  of  the  stress  series  is  plotted  in  Fig.  3.16  for  a  hardening 
exponent  of  n=10.  The  first  term  is  given  to  compare  the  influence  the  higher 
order  terms  have  on  the  results.  The  full  field  finite  element  results  for  each 
specimen  are  given  along  with  the  six  term  analytic  solution  in  Figs.  3.17-3.28 
for  a  hardening  exponent  of  n=10.  Tables  3.4-3.8  contain  the  general  power 


Table  3.1 —  Eigenvalue  ordering  for  plane  strain  and  a  hardening 
exponent  of  n=5. 


Eigenvalue 

Order 

Eigenvalue 

Si 

Difference 

ASi=Si-s, 

Particular 

Solutions 

Type  of  Differential 
Equation 

1 

-0.166670 

0.000000 

•  •• 

HOMOGENEOUS-PLASTIC 

2 

0.054558 

0.221220 

... 

HOMOGENEOUS-PLASTIC 

3 

0.275780 

0.442450 

1 

PARTICULAR-PLASTIC 

4 

0.340720 

0.507380 

HOMOGENEOUS-PLASTIC 

5 

0.497010 

0.663670 

2 

PARTICULAR-PLASTIC 

6 

0.500000 

0.666670 

1 

PARTICULAR-ELASTIC 

7 

0.561940 

0.728610 

1 

PARTICULAR-PLASTJC 

8 

0.718230 

0.884900 

4 

PARTICULAR-PLASTIC 

9 

0.721220 

0.887890 

1 

PARTICULAR-ELASTIC 

-100 

2.547700 

2.714400 

... 

• 

HOMOGENEOUS-PLASTIC 

Table  3.2 —  Eigenvalue  ordering  for  plane  strain  and  a  hardening 
exponent  of  n=10. 


Eigenvalue 

Order 

Eigenvalue 

Sj 

Difference 
As  j  =  Sj-s.| 

Particular 

Solutions 

Type  of  Differential 
Equation 

1 

-0.090909 

0.000000 

•  •• 

HOMOGENEOUS-PLASTIC 

2 

0.069766 

0.160680 

... 

HOMOGENEOUS-PLASTIC 

3 

0.230440 

0.321350 

1 

PARTICUUR-PLASTIC 

4 

0.269590 

0.360500 

... 

HOMOGENEOUS-PLASTIC 

5 

0.391120 

0.482030 

2 

PARTICULAR-PLASTIC 

6 

0.430270 

0.521180 

1 

PARTICULAR-PLASTIC 

7 

0.551790 

0.642700 

4 

PARTICULAR-PLASTIC 

8 

0.590950 

0.681850 

3 

PARTICULAR-PLASTIC 

9 

0.630100 

0.721010 

7 

PARTICULAR-PLASTIC 

10 

0.712470 

0.803380 

6 

PARTICULAR-PLASTIC 

11 

0.727270 

0.818180 

1 

PARTICULAR-ELASTIC 

-120 

2.065909 

1.975000 

... 

HOMOGENEOUS-PUSTIC 
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Table  3.3 —  Eigenvalue  ordering  for  plane  stress  and  a  hardening 
exponent  of  n=5. 


Eigenvalue 

Order 

Eigenvalue 

Sj 

Difference 

ASi=Si-s^ 

Particular 

Solutions 

Type  of  Differential 
Equation 

1 

-0.166670 

0.000000 

•  •• 

HOMOGENEOUS-PLASTIC 

2 

0.500000 

0.666670 

1 

PARTICULAR-ELASTIC 

3 

1.166700 

1.333300 

1 

PARTICUUR-ELASTIC 

4 

1.833300 

2.000000 

1 

PARTICULAR-ELASTIC 

5 

2.500000 

2.666700 

1 

PARTICUUR-ELASTIC 

6 

3.166700 

3.333300 

1 

PARTICUUR-EUSTIC 

7 

3.833300 

4.000000 

1 

PARTICUUR-EUSTIC 

8 

4.500000 

4.666700 

1 

PARTICUUR-EUSTIC 

9 

4.500000 

4.666700 

7 

PARTICUUR-PUSTIC 

10 

5.166700 

5.333300 

1 

PARTICUUR-EUSTIC 

FIG.  3.4—  Angular  stress  functions  for  n=10:  a)  field  1,  b)  field  2,  c)  field  3,  d)  field  4 
e)  field  5,  f)  field  6. 
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FIG.  3.5 —  Angular  strain  functions  for  n=10:  a)  field  1,  b)  field  2,  c)  field  3,  d)  field  4, 
e)  field  5,  f)  field  6. 
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FIG.  3.6 —  Angular  strain  functions  for  particular  solutions  to  differential  equations: 
a)  field  22,  b)  field  222,  c)  field  23,  d)  field  24. 


FIG.  3.10 —  Angular  stress  functions  for  n=5:  a)  field  1,  b)  field  2,  c)  field  3,  d)  field  4, 
e)  field  5,  f)  field  6. 
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FIG.  3.12 —  Angular  strain  functions  for  particular  solutions  to  differential  equations; 
a)  field  22,  b)  field  222,  c)  field  23,  d)  field  el. 


Angular  displacement  functions  for  n=5:  a)  field  1,  b)  field  2,  c)  field  3,  d)  field  4, 
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FIG.  3.14 —  Angular  displacement  functions  for  particular  solutions  to  differential  equations: 
a)  field  22,  b)  field  222,  c)  field  23,  d)  field  el. 
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15 —  Integration  constants  for  the  /  -  integral;  a)  field  1,  b)  field  2,  c)  field  3,  d)  field  4, 
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FIG.  3.17 —  Comparison  of  a  six  term  analytic  solution  with  finite 

element  predictions  for  the  in-plane  stresses  of  the  MBL: 

a)  Tja  =+0.3,  ct„,  =0.95(T  ,  b)  Tla  =-0.4,ct  =0.95ct  . 

'to  PL  0  '  I  0  PL  0 
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FIG.  3.18 —  Comparison  of  a  six  term  analytic  solution  with  finite 

element  predictions  for  the  in-plane  stresses  of  the  MBL: 

a)  77 <T  =-0.7,  =0.95(7  ,b)  Tie  =-0.7,  =0.85(7  . 

'  !  o  PL  0  '  I  0  PL  0 


FIG.  3.20 —  Comparison  of  a  six  term  analytic  solution  with  finite 
element  predictions  for  the  in-plane  stresses  of  the 
SENB  specimen  with  an  a/W=0.5:  a)  J /bcr^  =0.001462, 

b)  y/bcT^  =0.01111. 
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FIG.  3.22 —  Comparison  of  a  six  term  analytic  solution  with  finite 
element  predictions  for  the  in-plane  stresses  of  the 
CCT  specimen  with  an  a/W=0.1:  a)  J /bcr^  =0.000092, 
b)  7/5(7^=0.000421. 


FIG.  3.2S—  Comparison  of  a  six  term  analytic  solution  with  finite 
element  predictions  for  the  in-plane  stresses  of  the 
CCT  specimen  with  an  a/W^O-S:  a)  J /bor  =0.00071, 
b)  y/bcT^  =0.00259.  ° 
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FIG.  3.24 —  Comparison  of  a  six  term  analytic  solution  with  finite 
element  predictions  for  the  in-plane  stresses  of  the 
CCT  specimen  with  an  a/W=0.9:  a)  J  =0.001927, 

b)  Jb/CT^  =0.006356. 
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FIG.  3.25 —  Comparison  of  a  six  term  analytic  solution  with  finite 
element  predictions  for  the  in-plane  stresses  of  the 
DENT  specimen  with  an  a/W=0.1:  a)  J  =0.000319, 

b)  7/b(T^  =0.001039. 
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FIG.  3.26 —  Comparison  of  a  six  term  analytic  solution  with  finite 
element  predictions  for  the  in-plane  stresses  of  the 
DENT  specimen:  a)  a/W  =  0.5,  7/bcT^  =0.00245, 
b)  a/W  =  0.9, 7/b(T^  =0.001638. 
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FIG.  3.28 —  Comparison  of  a  six  term  analytic  solution  with  finite 
element  predictions  for  the  in-plane  stresses  of  the 
SENT  specimen:  a)  a/W  =  0.5,  7 /bcr^  =0.02230, 

b)  a  /  W  =  0.9,  J =0.024900. 
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series  coefficients  as  a  function  of  deformation  level  for  each  specimen 
geometry.  The  coefficients  were  obtained  as  a  function  of  position  within  an 
armulus  bounded  by  a  normalized  radius  of  1.0  <  r  <  15.  this  region  excludes 
the  finite  strain  zone  near  the  crack-tip.  It  was  found  that  at  low  to  medium 
deformations  levels,  the  value  of  coefficients  did  not  have  any  significant 
spatial  variation  within  this  annulus  except  for  very  near  the  crack  surfaces. 
The  coefficients  that  are  reported  were  determined  at  a  spatial  coordinate  of 
r  =  2.148  and  6  =  0.  Table  3.9  contains  the  coefficients,  deformation  levels  and 
location  for  the  analytic/numerical  comparison  given  in  Figs.  3.17-3.28. 

3.4  Influence  of  Second  and  Third  Williams  Terms 

To  insure  that  the  boundaries  in  the  MBL  had  no  influence  on  the  stress 
field  at  the  crack-tip,  the  displacements  associated  with  the  first  Williams  term 
were  applied  for  two  different  values  of  applied  K  .  A  comparison  of  the 
normalized  hoop  stress  are  given  in  Fig.  3.29.  The  stress  fields  after  being 
normalized  by  Oq  are  self  similar  in  nature  indicating  that  the  stresses  at  the 
CTack-tip  have  not  felt  the  influence  of  the  boundaries.  Further,  the  /-integral 
values  obtained  from  the  numerical  study  agree  to  within  a  one  half  of  one 
percent  to  that  of  the  applied  /  value. 

A  comparison  of  the  crack-tip  stresses  for  the  second  and  third 
Williams  terms  are  given  in  Figs.  3.30  and  3.31.  The  MBL  was  ran  first  with 
the  third  Williams  term  {S/Of,)  set  to  zero  while  varying  the  second  term 
(T/cTg)  and  then  the  second  term  was  set  to  zero  while  the  third  term  was 
varied.  The  values  of  the  second  and  third  coefficients  reported  in  these 
figures  correspond  to  plastic  zones  that  measure  approximately  8%  of  the 
MBL  radius.  The  first  fringe  plot  in  each  of  the  figures  shows  the  hoop  stress 
within  the  annulus  of  interest.  The  second  plot  in  the  figures  shows  the 
portion  of  the  hoop  stress  which  is  associated  with  the  second  and  third  term, 
respectfully.  The  last  fringe  plot  shows  the  global  hoop  stress  associated  with 
the  second  and  third  term,  respectfully.  The  contribution  of  the  second  and 
third  Williams  term  were  determined  by  subtracting  the  full  field  solutions  for 
T /(Tfl,  S/<jg  =  0  from  the  cases  for  T =  a,  S/a^  =  p. 
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TABLE  3.9 —  Coefficients  and  radial  value  used  to  compare  finite 
element  results  with  analytic  solution. 


Specimen 

Geometry 

Load 

r 

J 

K2 

K4 

T/o 

o 

(J/cm^  ) 

low  yield 

-0.70 

2.900 

1.32 

-0.0480 

-0.4750 

MBL 

high  yield 

-0.70 

2.900 

1.32 

-0.0480 

-0.4550 

high  yield 

-0.40 

2.867 

1.32 

-0.0280 

-0.4932 

high  yield 

0.30 

2.862 

1.32 

0.0200 

-0.1545 

(KN) 

J/ba 

0 

a/w  =  0. 1 

80.00 

2.286 

0.0005300 

-0.0230 

-0.4600 

142.22 

2.269 

0.0024778 

-0.0490 

-0.3350 

SENB 

a/w  =  0.5 

26.67 

2.108 

0.0014623 

0.0016 

-0.3763 

57.78 

2.174 

0.0111083 

-0.0034 

-0.5621 

a/w  =  0.9 

1.78 

2.101 

0.0038808 

0.0079 

-0.2617 

2.67 

2.075 

0.0320583 

0.0052 

-0.6193 

(MPa) 

a/w  =  0. 1 

155.17 

2.320 

0.0000920 

-0.0250 

-0.4920 

310.34 

2.149 

0.0004216 

-0.0510 

-0.3830 

CCT 

a/w  =  0.5 

120.69 

2.333 

0.0007063 

-0.0242 

-0.4924 

211.21 

2.044 

0.0025883 

-0.0496 

-0.3857 

a/w  =  0.9 

30.17 

2.312 

0.0019275 

-0.0260 

-0.4830 

47.41 

2.282 

0.0063567 

-0.0560 

-0.2880 

a/w  =  0.1 

256.03 

2.108 

0.0003194 

-0.0235 

-0.4471 

398.28 

2.054 

0.0010398 

-0.0485 

-0.3217 

DENT 

a/w  =  0.5 

224.14 

2.157 

0.0024533 

-0.0206 

-0.4570 

a/w  =  0.9 

34.48 

2.322 

0.0016383 

-0.0029 

-0.4112 

a/w  =  0.1 

241.38 

2.062 

0.0003265 

-0.0235 

-0.4773 

379.31 

2.037 

0.0012056 

-0.0503 

-0.3386 

SENT 

a/w  =  0.5 

168.97 

2.138 

0.0223000 

-0.0278 

-0.4695 

a/w  =  0.9 

4.17 

2.032 

0.0249000 

0.0029 

-0.5879 
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T/Oo  =  0.00 

S/ao  =  0.00  Kappiied  =  27.48  MPa  ml/2 


3.5472 

3.3067 

3.0662 

2.8257 

2.5852 

2.3447 

2.1042 

1.8637 

1.6232 

1.3827 

1.1422 

0.9017 

0.6612 

0.4207 

0.1802 

-0.0603 


Kapplied  =  54.96  MPa  ml/2 


r=10 


FIG.  3.29 —  Hoop  stress  for  first  term  of  modified  boundary  layer: 

a)  Kappiied  =  27.48  MPa  ml/2 1)  Kappiied  =  54.96  MPa  ml/2. 
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v/‘  ^/rit 
Kv.*  w 


<j09'<Jo  = -0-0599 
Ao0^ao=  0.2407 


-0.0109 

0.0009 


Oee/Oo  =  -0.0950 
Aae0/ao=  0.0334 


r  =  31648 


FIG.  3.30 —  Modified  boundary  layer  for  first  two  Williams  terms: 

a)  normalized  hoop  stress,  b)  contribution  from  second 
term  near  crack-tip,  c)  global  contribution. 
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=  -0.0950 
Agq^Oq  =  0.0334 


FIG.  3.31 —  Modified  boundary  layer  for  first  and  third  Williams  terms: 

a)  normalized  hoop  stress,  b)  contribution  from  third  term 
near  crack-tip,  c)  global  contribution. 
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3.5  Assumption  in  Constitutive  Equation  for  the  Out-of-Plane  Stress 

The  analytical  derivation  for  plane  strain  presumed  that  the  out-of¬ 
plane  stress  (Tjj  was  related  to  the  in  plane  stresses  by 

C^{r,e)^l{r,e)  [Grr{r,e)-i-C7Qg{r,e)] 


where 


X(r,0)  = 


v+|or'('-.e) 

[l-i-aar^(r,e)]  ’ 


The  assumption  was  made  that  within  a  fully  developed  plastic  region,  as 
would  exist  at  the  crack-tip,  v  «  resulting  in  a  limiting  value  of  X  0.5, 
Figs.  3.32  and  3.33  show  the  value  of  X(r,  6)  obtained  from  the  MBL 
results  for  different  values  of  the  second  and  third  Williams  term.  The 
annulus  of  this  study  is  contained  within  a  normalized  radius  of  ten.  The 
results  indicate  that  the  region  in  which  X  is  approaching  the  limiting  value  of 
0.5  depends  on  the  second  Williams  term.  Setting  the  second  Williams  term 
equal  to  zero  results  in  the  second  fringe  plot  in  Fig.  3.32.  The  region  in  which 
X->0.5  is  at  Applying  negative  T/Cg  values  rotates  this  region 

towards  the  aack  plane  (first  fringe  plot)  while  positive  values  of  T/Og  rotate 
the  region  towards  the  crack  faces  (third  fringe  plot).  The  region  on  the  crack 
plane  that  corresponds  to  values  of  X  >  0.4  resides  within  a  normalized  radius 
of  approximately  four.  As  noted  previously,  the  third  Williams  term  has  little 
influence  on  the  crack-tip  region  for  values  of  S/Og  which  do  not  cause  plastic 
flow  beyond  10%  of  the  MBL  radius.  Fig.  3.33  shows  X(r,&)  for  values  of 
S/(7g  that  satisfy  the  above  condition  in  all  three  fringe  plots  the  region  is  at 

e^n/2. 


3.6  Evolution  of  the  Second  and  Third  Independent  Amplitudes 


The  evolution  of  the  second  (Kz  )  and  third  (K4  )  independent 
parameters  (Tables  3.4-3,8)  are  shown  in  Figs.  3.34-3.36  as  a  function  of  the 
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FIG.  3.33 —  Lambda  as  a  function  of  S/ Oq  and  position  for  the  MBL, 
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deformation  level  (//bOo)  and  a/W  for  the  SENB,  CCT,  DENT,  and  SENT 
specimens,  respectfully.  These  data  are  compared  with  the  Williams  two- 
parameter  surface  obtained  for  the  MBL.  Dashed  lines  indicate  when  the 
finite  geometry  results  have  significantly  deviated  from  the  MBL  surface.  The 
two-parameter  surface  is  constructed  from  the  K2  and  K4  relation  obtained  for 
the  MBL.  The  Kz  and  K4  relation  can  be  extruded  in  the  vertical  direction  as  a 
result  of  the  self  similarity  in  the  normalized  stress  field  for  increasing 
deformation  level  as  shown  in  Fig.  3.29.  This  surface  by  the  definition  of  the 
boimdary  displacements  applied  to  the  MBL  is  orUy  dependent  on  two- 
parameters  {K,  T). 

The  evolution  of  the  second  independent  coefficient  for  each  of  the 
geometries  is  plotted  as  a  function  of  deformation  level  (J/hOo)  in  Figs.  3.37- 
3.39  for  a/W=0.1,  0.5,  and  0.9,  respectfully.  Also  included  is  the  value  of  the 
second  coefficient  obtained  for  the  SSY  model  (MBL  ->  T  =  0).  Fig.  3,40 
compares  the  relationship  between  the  second  and  third  independent 
parameters  for  the  second  and  then  the  third  Williams  term  MBL  results.  Two 
sets  of  the  second  Williams  term  MBL  are  shown.  The  first  corresponds  to  a 
linear,  elliptic,  power-law  constitutive  relation  with  a  proportional  limit  of 
85%  of  the  flow  stress.  The  second  of  the  MBL  results  is  based  on  a 
proportional  limit  of  95%  of  the  flow  stress.  For  values  of  the  proportional 
limit  less  than  about  92%  boimdary  plasticity  occurred  at  highly  negative 
values  of  T/Cg  Hence,  the  plastic  zone  was  influenced  by  the  proximity  of  the 
MBL  boundary.  At  values  greater  than  92%  the  plasticity  was  contained 
within  a  radius  that  was  10%  of  the  MBL  radius.  The  relationship  between  Kz 
and  K4  did  not  vary  for  highly  negative  values  of  T/Cg  values. 

The  relationship  between  Kz  and  K4  are  given  in  Figs.  3.41-3.44  for  each 
of  the  specimen  geometries.  The  two-parameter  MBL  results  based  on  the  T 
-stress  are  also  plotted.  Figs.  3.45-3.47  compare  the  second  and  third 
independent  parameter  relation  for  different  specimens  at  a/W  ratios  of  0.1, 
0.5,  and  0.9,  respectfully. 

The  magnitude  of  the  difference  field  as  given  in  Eq.  (1.10)  is 
determined  by  computing  the  difference  between  the  third  independent 
coefficient  K4  for  the  finite  geometry  specimens  with  that  of  the  second 
Williams  term  two-parameter  MBL  value  of  K 4  for  a  fixed  Kz  value.  The 
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FIG.  3.39 —  Comparison  of  the  second  independent  amplitude  with 
that  obtained  for  the  SSY  model  for  the  specimei\s  with 
a/W  =0.9. 
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FIG.  3.40 —  Comparison  of  the  normalized  amplitudes  for  the  MBL 
between  the  first  and  second  Williams  terms  and 
different  proportional  limits  in  the  constitutive  relation. 
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FIG.  3.43 —  Comparison  of  the  normalized  amplitudes  for  the  DENT 
specimen  with  those  for  the  two-parameter  MBL. 


45 —  Comparison  of  the  normalized  amplitudes  for  the 
specimens  having  a/W  =  0.1  with  those  for  the  two- 
parameter  MBL. 


Comparison  of  the  normalized  amplitudes  for  the 
specimens  having  a/W  =  0.5  with  those  for  the  two- 
parameter  MBL. 


Comparison  of  the  normalized  amplitudes  for  the 
specimens  having  a/W  =  0.9  with  those  for  the  two- 
parameter  MBL. 


FIG.  3.49 —  Difference  in  the  third  independent  coefficients  between 
the  CCT  and  MBL. 


FIG.  3.52 —  Difference  between  the  third  independent  coefficients  for 
the  specimens  having  a/W  =  0.1  with  those  for  the  two- 
parameter  MBL. 


FIG.  3.54 —  Difference  between  the  third  independent  coefficients  for 
the  specimens  having  a/W  =  0.9  with  those  for  the  two- 
parameter  MBL. 
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percent  error  is  given  in  Figs.  3.48-3.51  for  the  SEND,  CCT,  DENT,  and  SENT, 
respectfully  Figs.  3.52-3.54  show  the  percent  error  between  K4  for  a/W  values 
of  0.1,  0.5  and  0.9,  respectfully.  A  dashed  line  has  been  added  to  indicate 
when  the  finite  body  stresses  have  deviated  more  than  20%  from  the  two- 
parameter  MBL  stresses.  The  corresponding  deformation  levels  as  measured 
by  J/hCo  are  tabulated  in  Table  3.10. 

The  error  in  the  stress  field  as  given  in  Eq.  (1.10)  associated  with  a  20% 
difference  between  the  finite  geometry  third  parameter  and  the  MBL  third 
parameter  are  given  in  Figs.  3.55  and  3.56  for  the  SENB  and  CCT  specimens 
having  a/W=0.5.  Fig.  3.55  shows  the  corresponding  errors  along  the  crack 
plane  while  Fig.  3.56  shows  the  error  in  the  angular  direction  at  a  normalized 
radius  of  five. 
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TABLE  3.10 —  Deformation  levels  for  20%  difference  between  the  stress 
fields  in  the  two-parameter  MBL  and  the  finite  geometries. 


J/bq> 

a/w 

0.1 

0.5 

0.9 

SENS 

0.00198 

0.00026 

0.00772 

CCT 

0.00052 

0.00259 

0.00457 

DENT 

0.00061 

0.00195 

— 

SENT 

0.00079 

— 

0.00273 

FIG.  3.55 —  Error  field  as  a  function  of  the  radial  coordinate  for  the 
SENB  and  CCT  specimens  at  deformations 
corresponding  to  a  20%  error  in  the  third  independent 
parameter  JQ. 


6  (radians) 


FIG.  3.56 —  Error  field  as  a  function  of  theta  for  the  SENB  and  CCT 
specimens  at  deformations  corresponding  to  a  20%  error  in 
the  third  independent  parameter  K4. 
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CHAPTER  IV 
DISCUSSION 


In  this  work  we  determine  when  the  stresses  in  the  plastic  zone  deviate 
from  those  characterized  by  a  two-parameter  methodology.  In  so  doing,  we 
utilize  the  general  power  series  representation  for  the  stress  potential.  As  far 
as  the  analytic  analysis  for  the  series  representation  of  the  crack-tip  stress 
fields,  it  is  impossible  to  transition  from  the  plastic  region  to  the  elastic  region 
without  assuming  the  material  is  elastic  incompressible  (i.e., 
v  =  0.3  ->X(r,0)  =  O.5).  The  sixth  field  for  the  case  of  plane  strain  with  a 
hardening  coefficient  of  n=5  is  the  first  elastic  term  (Table  3.1).  Fig.  3.15f 
shows  the  sixth  integration  constant  (Ig-  Appendix  F)  for  the  /-integral.  The 

constant  must  integrate  to  zero  to  preserve  the  path  independence  of  /.  The 
integration  constant  is  shown  for  both  v=0.3  and  v=0.5.  The  integral  is  not 
zero  for  v=0.3  leaving  the  /-integral  path  dependent.  This  discrepancy  arises 
as  a  result  of  setting  X(r,0)  =  0.5  in  the  third  constitutive  equation,  Eq,  (2.51). 
This  assumption  was  necessary  to  make  the  problem  mathematically  tractable. 
The  validity  of  this  assumption  was  addressed  by  looking  at  the  distribution 
of  X(r,  6)  in  the  MBL.  Figs.  3.32  and  3.33  indicate  the  region  were  X{r,  6) -¥0.5 
is  at  6  =it/2  and  within  a  normalized  radius  of  approximately  four.  Hence, 
some  error  is  introduced  with  this  assumption. 

The  comparison  of  the  analytic  solutions  with  the  finite  element  results 
in  Figs.  3.17-3.28  indicate  that  the  six  term  series  representation  for  the  stress 
field  seems  to  do  an  adequate  job  in  describing  the  stress  field  in  the  annular 
region  bounded  by  a  normalized  radius  of  five. 

The  two-parameter  solution  that  we  base  our  analysis  is  that  for  the 
MBL  with  a  two  term  Williams  displacement  field  applied  on  its  elastic 
boundary.  Comparison  of  the  second  fringe  plot  in  Figs.  3.30  and  3.31  clearly 
indicates  that  the  second  Williams  term  has  a  greater  influence  on  the  stresses 
in  the  plastic  zone  then  that  of  the  third  term.  The  reason  for  this  is  that  the 
stress  associated  with  the  second  Williams  term  has  no  dependence  on  the 
radial  coordinate  while  the  third  term  has  a  dependence.  Hence,  the 


137 


nearer  the  crack-tip,  the  smaller  the  stress  associated  with  the  third  term.  In 
general,  this  would  indicate  that  if  the  elastic  field  outside  the  plastic  region 
can  be  characterized  by  the  Williams  solution,  then  the  higher  order  terms 
have  a  diminishing  effect  on  the  crack-tip  plasticity. 

The  evolution  of  the  independent  coefficients  for  the  series 
representation  of  the  stresses  in  Figs.  3.34-3.36  dearly  indicate  that  the  shallow 
cracked  geometries  (a/W=0.1)  maintain  a  two-parameter  field  description 
longer  then  that  of  some  of  the  medium  and  deeply  cracked  geometries.  The 
second  parameter  can  be  incorporated  in  a  failure  locus  at  the  onset  of 
deavage  by  determining  its  value  at  Jc  in  Figs.  3.37-3.39. 

The  relation  between  the  second  and  third  independent  coefficients  in 
the  analytic  elastic-plastic  solution  for  both  the  MBL  and  finite  geometries  are 
shown  in  Figs.  3.45-3.47.  For  a/W=0.1,  the  CCT  geometry  maintains  a  two- 
parameter  description  to  higher  deformation  levels  then  those  assodated  with 
the  SENB,  DENT,  and  SENT  specimens.  This  result  can  best  be  seen  in  Figs. 
3.34  and  3.45  were  the  coeffident  for  the  difference  field  between  the  MBL  and 
finite  geometries  with  an  a/W  =  0.1  are  shown. 

The  CCT  also  maintains  a  two-parameter  field  to  higher  deformation 
levels  then  the  other  finite  geometries  for  both  a/W  =  0.5  and  0.9.  The  SENB 
specimen,  as  shown  in  Figs.  3.35,  3.47,  and  3.40,  deviates  from  the  two- 
parameter  description  at  quite  low  deformation  levels  when  a/W  =  0.5  and 
0.9,  respectfully.  The  SENT  behaves  similar  in  nature  to  the  SENB  specimens 
at  a/W  =  0.9  because  of  the  presence  of  the  large  global  bending  field  common 
in  the  SENB  specimens  at  high  a/W  values.  One  difference  between  the  SENB 
specimen  with  an  a/W  =  0.9  and  the  remaining  geometries,  is  that  there  is  an 
increase  in  constraint  before  the  global  bending  field  becomes  dominate.  The 
bending  field  can  be  seen  to  dominate  in  the  finite  element  results  shown  in 
Figs.  3.20  and  3.21  for  a/W  =  0.5  and  a/W  =  0.9,  respectfully. 

The  third  Williams  term  was  originally  added  to  the  boundary  of  the 
MBL  in  an  attempt  to  account  for  the  global  bending  field  present  in  the  SENB 
(a/W  =  0.5  and  0.9)  and  SENT  (a/W  =  0.9)  finite  geometries.  The  third  term  as 
previously  mentioned  had  little  influence  on  the  stress  field  within  the  plastic 
zone  for  S/ Oq  values  for  which  the  plastic  zone  did  not  feel  the  influence  of 
the  bovmdaries  of  the  MBL. 
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The  results  for  the  difference  fields  in  Figs.  3.48-3.54  show  that  of  the 
finite  geometries  modeled,  all  but  the  SENB  (a/W  =  0.5, 0.9)  and  SENT  (a/W  = 
0.9)  have  stress  field  within  the  plastic  zone  that  can  be  characterized  by  two 
parameters  at  low  to  moderate  deformation  levels.  What  are  the  implications 
drawn  from-  these  results?  The  best  specimens  for  two-parameter 
characterization  of  the  stress  field  would  appear  to  be  the  CCT,  SENT  (a/W  = 
0.1  and  0.5)  and  CCT  (a/W  =  0.9).  Although,  since  the  deep  notched  SENB 
and  SENT  specimens  exhibit  single-parameter  dominance  by  remaining  closer 
to  the  SSY  solution  at  low  to  moderate  deformation  levels,  the  difference  in  a 
higher  parameter  description  may  be  insignificant.  The  error  in  the  difference 
field  (Eq.  1.10)  for  20%  deviation  from  the  SSY  surface  for  the  SENB  (a/W=0.5) 
is  less  then  that  for  the  corresponding  CCT  specimen  as  seen  in  Figs.  3.55  and 
3.56. 
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CHAPTER  V 

SUMMARY  AND  CONCLUSIONS 


In  the  course  of  this  study,  the  fifth  and  sixth  angular  functions  and 
corresponding  radial  exponents  for  the  generalized  power  series 
representation  of  the  stress  potential  were  determined  for  a  Ramberg-Osgood 
material  having  hardening  exponents  of  n=5  and  10.  The  leading  coefficients 
for  each  term  were  determined  for  four  common  laboratory  specimen 
geometries  as  a  function  of  deformation  level  by  matching  the  analytic 
solution  with  full  field  finite  element  results. 

The  general  power  series  representation  of  the  stress  field  can  only  be 
determined  up  to  the  fifth  term  for  a  hardening  exponent  of  n=5  and  up  to  the 
tenth  term  for  a  hardening  exponent  of  n=10.  The  next  terms  in  the  series 
correspond  to  the  elastic  portion  of  the  strain.  The  angular  coefficient  for  the 
/-integral  does  not  integrate  to  zero  for  these  terms  unless  an  elastic 
incompressible  material  is  assumed.  This  difficulty  arises  due  to  the  fact  that 
the  effective  plastic  poisson  ratio  Ain  the  third  constitutive  equation  is 
assumed  to  be  0.5  for  the  case  of  plane  strain.  When  an  elastic  incompressible 
material  having  a  poisson  ratio  of  0.5  is  assumed,  the  sixth  angular  coefficient 
integrates  to  zero  preserving  the  path  independence  of  the  /-integral. 

A  criterion  was  established  by  which  it  is  possible  to  determine  when 
the  near  tip  stress  fields  deviate  from  characterization  by  two-parameter 
descriptions.  Comparisons  of  the  evolution  of  the  independent  parameters  in 
the  general  power  series  representation  of  the  stress  field  for  the  laboratory 
specimen  geometries  with  those  for  the  two-parameter  MBL  solution  indicate 
that  all  but  the  SENS  (a/W  =  0.5  and  0.9)  and  SENT  (a/W  =  0.9)  lend 
themselves  to  a  two-parameter  description  for  low  to  moderate  deformation 
levels.  The  SENB  and  SENT  geometries  that  do  not  stay  in  the  two-parameter 
surface  of  the  MBL  at  low  and  moderate  deformation  levels  however  stay  near 
the  SSY  solution.  Consequently,  a  single  parameter  description  is  a  good 
approximation.  As  the  a/W  ratio  of  the  specimens  increases,  we  find  the  path 
determined  by  the  first  three  independent  coefficients  leaves  the  two- 
parameter  MBL  surface  at  lower  deformation  levels. 
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The  results  indicate  that  for  a  20%  difference  in  the  third  independent 
coefficient  for  the  SENB  and  CCT  as  compared  with  the  MBL  geometry,  the 
error  in  the  stress  field  as  predicted  by  the  general  power  series  representation 
of  the  stress  field  is  less  for  the  SENB  as  compared  to  the  CCT  specimen.  The 
reason  for  this  is  that  the  SENB  specimen  is  at  a  lower  deformation  level  then 
the  CCT  specimen  when  there  is  a  20%  difference  between  the  two-parameter 
surface  and  the  finite  geometry  results. 

The  laboratory  specimens  studied  in  this  work  did  not  include  the 
compact  tension  specimen  (CT).  Future  work  along  the  same  line  as  the 
present  research  could  include  determining  the  evolution  path  of  the  first 
three  independent  coefficients  in  the  general  power  series  representation  for 
the  stress  field.  Because  this  specimen  demonstrates  predominately  a  global 
bending  field  on  the  crack  plane  it  is  fortuitous  to  expect  the  results  to  be 
similar  to  that  of  the  SENB  specimen. 
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Q  ************************************************************** 

SUBROUTINE  UMAT (STRESS, STATEV, DDSDDE, SSE, SPD, SCD, 

Sc  RPL ,  DDSDDT ,  DRPLDE ,  DRPLDT , 

&  STRAN, DSTRAN, TIME, DTIME, TEMP, DTEMP, PREDEF , DPRED , CMNAME, 
t  NDI , NSHR, NTENS , NSTATV, PROPS, NPROPS ,  COORDS , DROT, PNEWDT, 

&  CELENT) 

IMPLICIT  REAL*8  (A-H,0-Z) 

PARAMETER {NPRECD=2 ) 

CHARACTER* 8  CMNAME 

DIMENSION  STRESS (NTENS ) , STATEV (NSTATV) , 

£c  DDSDDE  ( NTENS ,  NTENS )  ,  DDSDDT  ( NTENS )  ,  DRPLDE  ( NTENS )  , 

&  STRAN ( NTENS ) ,DSTRAN( NTENS )  ,TIME(2) , PREDEF (1) ,DPRED(1) , 

&  PROPS (NPROPS ) , COORDS ( 3 ) ,  DROT (3,3) 

DIMENSION  AA(4,4) ,BB(4,4) ,FF(4) ,DD(4) 

DIMENSION  STRANEW ( 6 ) , EDEV ( 6 ) 

Define  material  constants 

SO=PROPS(l) 

EO=PROPS(2) 

RN=PROPS(3) 

RNU=PROPS(4) 

EMOD=SO/EO 

Constants  for  the  transition  region 

ONE=l. 000000000000000000000000 
ZERO=0 . 00000000000000000000000 
TOLERANCE=0 . 0000000001 
PI=4.0*ATAN(ONE) 

PI2=PI/2.0 
RK2=1.1 
RK1=0.85 

EE_1=2 . 0* ( 1+RNU) *RKl*EO/3 . 0 
EE_2 = ( RK2 * *RN+ 2 . 0 * ( 1 . 0+RNU) *RK2 / 3 . 0 ) *EO 
C1=RK1 
C2=1.0 

C3=RK2+RK2**RN 
C4=l . 0+RN*RK2** (RN-1 .0) 

QNC= (C2/C4*RK2-RK1) / (C2/C4-1 . 0) -0 . 01 
R2=( (RK2-QNC)**2.0*(1.0-C2**2.0/C4**2.0) 
t  / (1.0-C2**2.0/C4**2.0*(RK2-QNC)**2.0 
&  /(RK1-QNC)**2.0))**0.5 

C 

III=1 

c 

100  FF(1)=-(R2**2.0*( (RK1-QNC)**2.0-C2**2.0/C4**2.0*(RK2QNC)**2.0) 

&  -(RK1-QNC)**2.0*(RK2-QNC)**2.0*(1.0-C2**2.0/C4**2.0) ) 

FF(2)=-(C4*(R2**2.0-(RK2-QNC)**2.0)- 
&  C4*(R2**2.0-(RK2-QNC)**2.0)**0.5* 

fit  (R2**2.0-(RK1-QNC)**2.0)**0.5-(C1-C3)* (RK2-QNC) ) 
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AA(1, 1)=-2.0*R2**2.0*(RK1-QNC) 

Sc  +2.0*R2**2.0*C2**2.0/C4**2.0*{RK2-QNC) 

Sc  +2.0*(1.0-C2**2.0/C4**2.0)  *(RK1-QNC)*{RK2-QNC) 

Sc  *(  (RK1-QNC)  +  {RK2-QNC) ) 

AA{1,2)=(RK1-QNC)**2.0-C2**2.0/C4**2.0*{RK2-QNC)**2.0 
AA{2,1)=2.0*C4*(RK2-QNC)-C4* 

Sc  (  (RK1-QNC)*(R2**2.0-(RK2-QNC)  **2.0)  **0,5 

Sc  '  /(R2**2.0-(RK1-QNC)**2.0)**0.5 

Sc  +  (RK2-QNC)  *  (R2**2 . 0-  (RKl-QNC)  **2 . 0)  **0 . 5 

Sc  /(R2**2.0-(RK2-QNC)**2.0)**0.5)  +  (C1-C3) 

AA{2,2)=:2.0*C4*R2-C4*R2* 

Sc  ( (R2**2.0-(RK2-QNC)**2.0)**0.5 

Sc  /{R2**2.0-{RK1-QNC)**2.0)**0.5 

Sc  +{R2**2.0- (RKl-QNC)  **2.0)  **0.5 

Sc  /{R2**2.0-(RK2-QNC)**2.0)**0.5) 

C 

BB(1,1)=AA(2,2)/{AA{1,1)*AA{2,2)-AA(1,2)*AA{2,1) ) 
BB(1,2)=-AA(1,2)/(AA(1.1)*AA{2.2)-AA(1,2)*AA(2,1) ) 
BB(2,1)=-AA(2,1)/(AA(1,1)*AA(2,2)-AA{1,2)*AA(2,1) ) 
BB(2,2)=AA(1,1)/(AA(1,1)*AA(2,2)-AA{1,2)*AA(2,1) ) 

C 

DD(1)=BB(1,1)*FF(1)+BB(1,2)*FF(2) 

DD(2)=BB(2, 1)*FF(1)+BB(2,2)*FF(2) 

C 

RMAXDIFF=ZERO 

IF(ABS(DD(1) ) .GT,RMAXDIFF)RMAXDIFF=ABS(DD(1) ) 
IF{ABS(DD(2) ) .GT.RMAXDIFF)RMAXDIFF=ABS(DD(2) ) 

C 

QNC2aQNC+DD(l) 

R22=R2+DD(2) 

C 

QNC=QNC2 

R2=R22 

C 

IF ( RMAXDIFF . GT . TOLERANCE ) THEN 
IF(III.GT.300)THEN 

WRITE (7,*)  'Number  of  ellipse*, 

S;  '  iterations  to  high!  ' 

STOP 
ENDIF 
III=III+1 
GOTO  100 
ENDIF 
C 

R2=( (RK2-QNC) **2.0* (1 . 0-C2**2 . 0/C4**2 . 0 ) 

Sc  /(1.0-C2**2.0/C4**2.0*(RK2-QNC)**2.0 
Sc  /{RK1-QNC)**2.0))**0.5 
Rl=C4*R2*(1.0-{ {RK2-QNC)/R2)**2.0)**0.5 
Sc  /  {  {RK2-QNC) /R2) 

EENC=C1+R1* { 1 . 0- ( (RKl-QNC) /R2) **2 . 0 ) **0 . 5 
C 

C  Calculate  total  strain  from  begining  strain  and  strain 

C  increment 

C 


DO  I=1,NTENS,1 


nnnnnnnnnnn  non  n  n  n  n  ooo  oon 
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STRANEW ( I ) =STRAN ( I ) +DSTRAN ( I ) 

ENDDO 

Calculate  volumetric  strain  and  deviatoric  strain 

EV= ( STRANEW ( 1 ) +STRANEW ( 2 ) +STRANEW ( 3 ) ) /3 . 0 
DO  I=1,ND1,1 

EDEV ( I ) =STRANEW ( I ) -EV 
ENDDO 

DO  I=1,NSHR,1 

EDEV(NDI+I) =STRANEW{NDI+I) 

ENDDO 

Calculate  equivalent  strain 

EE=0. 00000 
DO  I=1,ND1,1 

EE=EE+EDEV(I)**2.0 

ENDDO 

DO  I=1,NSHR,1 

EE=EE+2 .0* (EDEV(NDI+I) /2 . 0) **2 . 0 
ENDDO 

EE={2.0*EE/3,0)**(0,5) 

Zero  stress  array  (given  strains  there  is  only  one  stress 
tensor) 

DO  I=1,NTENS,1 
STRESS ( I ) =2ERO 
ENDDO 

Zero  jacobian 

DO  I=1,NTENS,1 
DO  J=1,NTENS,1 

DDSDDEd,  J)=ZERO 
ENDDO 
ENDDO 

Linear  elastic  segment 

****ir****************lr****itiriririr*it***********lt*ir**ir1rir***1tirir**** 

IF(EE.LE.EE_1)THEN 

Equivalent  stress 

SENEW=3 . 0/2 . 0*EMOD/ ( 1 . 0+RNU) *EE 

Stress  tensor 

DO  I=1,NDI,1 
DO  J=1,NDI,1 

STRESS { I ) =STRESS ( I ) +DDSDDE ( I ,  J ) *  STRANEW { J ) 

ENDDO 


onoooono  oooo  non  n  n  n 
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ENDDO 

DO  I=1,NSHR,1 

STRESS ( I+NDI ) =DDSDDE ( I+NDI , I+NDI ) *STRANEW ( I+NDI ) 
ENDDO 

Elastic  jacobian 

DO  I=1,NDI,1 
DO  J=1,NDI,1 

DDSDDEd,  J)  =EMOD/ (1.0+RNU)  *RNU/  ( 1 . 0-2 . 0*RNU) 
ENDDO 
ENDDO 

DO  I=1,NDI,1 

DDSDDE (1,1) =DDSDDE (1,1) +EMOD/ ( 1 . 0+RNU ) 

ENDDO 

DO  I=1,NSHR,1 

DDSDDE ( I+NDl ,  I+NDI ) =EMOD/ ( 1 . 0+RNU) /2 . 0 
ENDDO 


Update  state  variables 

STATEV ( 1 ) =EE 
STATEV(2)=SENEW 
STATEV ( 3 ) =EE*SENEW 

Update  strain  energy,  plastic  dissipation,  and  creep 
dissipation 

PRESS=2ERO 
DO  I=1,NDI,1 

PRESS=PRESS+STRESS ( I ) /3 . 0 
ENDDO 

SSE= ( (1+RNU) *SENEW**2 .0/3 . 0+3 . 0* ( 1 .0-2 . 0*RNU) 

&  *PRESS**2.0/2.0)/EMOD 

SPD=ZERO 
SCD=ZERO 

Transitional  Segment 

*******************************+****************************** 


ELSEIF( (EE.GT.EE_1) .AND. (EE.LT.EE_2) )THEN 


Equivalent  stress  by  Newtons  method 


III=1 

SE=STATEV(2) 

SESO=SE/SO 
IF ( SESO . LT . RKl ) THEN 
SES0=RK1 
ENDIF 
EEEO=EE/EO 

D=EENC-R1 * ( 1 . 0 - ( ( SESO-QNC )/R2)**2.0)**0.5+ 
&  (2.0*RNU-1.0)/3.0*SESO-EEEO 


o  o  o  o  o  o 
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200  STIFF=R1/R2* (SESO-QNC) /R2/ (1 . 0- { (SESO-QNC) /R2 ) **2 . 0) **0 . 5 
&  +(2.0*RNU-1.0)/3.0 

SESONEW=SESO-D/ STIFF 
IF ( SESONEW . GT . RK2 ) SESONEW=RK2 

DNEW=EENC-R1* (1 . 0- ( (SESONEW-QNC) /R2) **2 . 0) **0 , 5+ 

&  {2.0*RNU-1.0)/3.0*SESONEW-EEEO 

IF(III.GT.300)THEN 

WRITE(7,*)  ‘To  many  iterations  in  the  material  subroutine' 
STOP 
ENDIF 

IF(ABS(DNEW) .LT. TOLERANCE) THEN 
SENEW=SESONEW* SO 
ELSE 

IF ( SESONEW . LT . RKl ) THEN 
SESONEW= (RKl+SESO) /2 . 0 

DNEW=EENC-R1 * ( 1 . 0 - ( { SESONEW-QNC )/R2)**2.0)**0.5+ 

&  (2.0*RNU-1.0)/3.0*SESONEW-EEEO 

ENDIF 

SESO=SESONEW 
D=DNEW 
III=III+1 
GOTO  200 
ENDIF 

Stress  tensor 

IF ( EE. NE. ZERO) THEN 
DO  I=1,NDI,1 

STRESS ( I ) =2 . 0/3 . 0*SENEW/EE*EDEV ( I ) 

&  +EMOD/{1.0-2.0*RNU)*EV 

ENDDO 

DO  I=1,NSHR,1 

STRESS ( I+NDI ) =2 . 0/3 . 0*SENEW/EE*EDEV ( I+NDI ) /2 . 0 
ENDDO 
ENDIF 

Transitional  jacobian 

STIFF=R1/R2* (SESONEW-QNC) /R2/ 

Sc  (1.0-{ (SESONEW-QNC) /R2) **2.0) **0.5+ (2. 0*RNU-1.0)/3.0 
F=2 . 0/3 . 0*SENEW/EE 

G=1.0/3.0*(EMOD/(1.0-2.0*RNU)-2.0/3.0*SENEW/EE) 

IP (EE.NE. ZERO) THEN 

H=4 .0/9 . 0*SENEW/EE**2 . 0* (EMOD/ (SENEW*STIFF) -1 . 0/EE) 

ELSE 
H=ZERO 
ENDIF 
C 

C  Calculate  first  part  of  jacobian  while  also  changing 

C  from  engineering  to  tensor  strains 

C 

DO  I=1,NDI+NSHR,1 
DO  J=1,NDI+NSHR, 1 

IF( (I.LE.NDI) .AND. (J.LE.NDI) )THEN 
DDSDDE ( I , J ) =H*EDEV ( I ) *EDEV ( J ) 


ooo  ooo  ono 
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ELSEIF( (I.LE.NDI) .AND. (J.GT.NDI) ) THEN 
DDSDDE ( I , J ) =H*EDEV ( I ) *EDEV { J) /2 . 0 
ELSEIF( (I.GT.NDI) .AND. (J.LE.NDI) )THEN 
DDSDDE ( I , J ) =H*EDEV ( I ) *EDEV ( J ) /2 . 0 
ELSEIF{ (I.GT.NDI) .AND. (J.GT.NDI) )THEN 
DDSDDEd,  J)=H*EDEV(I)*EDEV(J)/2.0/2.0 
ELSE 

WRITE (7, *) 'Did  not  add  first  jacobian' 

ENDIF 
ENDDO 
ENDDO 

Calculate  second  part  of  jacobian 

DO  I=1,NDI,1 
DO  J=1,NDI,1 

DDSDDE ( I , J ) =DDSDDE ( I , J ) +G 
ENDDO 
ENDDO 

Calculate  third  part  of  jacobian 

DO  I=1,NDI+NSHR,1 

DDSDDE (1,1) =DDSDDE ( I ,  I ) +F 
ENDDO 

Change  jacobian  to  engineering  strain  (2eps  =  gamma) 

DO  I=1,NDI+NSHR,1 
DO  J=1,NDI+NSHR,  1 

IF( (I.GT.NDI) .OR. (J.GT.NDI) ) THEN 
DDSDDEd,  J)=DDSDDE (I.  J)  /2.0 
ENDIF 
ENDDO 
ENDDO 

Update  state  variables 

STATEV(1)=EE 
STATEV(2)=SENEW 
STATEV ( 3 ) =EE*SENEW 

C  Update  strain  energy,  plastic  dissipation,  and  creep 

C  dissipation 

C 

PRESS=ZERO 
DO  I=1,NDI,1 

PRESS=PRESS+STRESS ( I ) /3 . 0 
ENDDO 

SS_1=RK1*S0 

THETA1=ASIN( (SS_l/SO-QNC) /R2) 

THETA2=ASIN( (SENEW/SO-QNC) /R2) 

SSE_PRESS=3 . 0* (1 . 0-2 . 0*RNU) /2 . 0*PRESS**2 . 0 
SSE_ELASTIC= ( 1 . 0+RNU) /3 . 0*SS_1**2 . 0 
SSE_TRANS= (Rl*R2*SO**2 . 0* (1 .0/2 . 0*THETA2-1 .0/4.0 
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&  *SIN{2.0*THETA2) 

fit  -QNC/R2*{1.0-{ (SENEW/S0-QNC)/R2)**2.0)**0.5) ) 

&  -(R1*R2*SO**2.0*(1. 0/2. 0*THETA1-1. 0/4.0 

&  *SIN(2.0*THETA1) 

Sc  -QNC/R2*(1.0-{ (SS_1/SO-QNC)/R2)**2.0)**0.5) ) 

SSE= (SSE_PRESS+SSE_ELASTIC+SSE_TRANS) /EMOD 
SPD=ZERO 
SCD=ZERO 
C 

Q  ★★★★♦**★**★*★**★**★★***★♦★★★**♦*******★★**★***★*★★★*★★*★***★*★ 

C  Power  law  Segment 

c 

ELSEIF ( EE . GE . EE_2 ) THEN 
C 

C  Equivalent  stress  by  Newtons  method 

C 

III=1 

TOLERANCE=0 . 0000000001 
SE=STATEV{2) 

SESO=SE/SO 

EEEO=EE/EO 

D=SESO**RN+2 . 0* ( 1+RNU) /3 . 0*SESO-EEEO 
300  ST1FF=RN*SES0** (RN-1 ,0)+2 . 0/3 . 0* (1 . 0+RNU) 

SESONEW=SESO-D/STIFF 

DNEW=SESONEW**RN+2 . 0* (1+RNU) /3 . 0*SESONEW-EEEO 
IF(III.GT,300)THEN 

WRITE(7,*)  'To  many  iterations  in  the  material  subroutine' 
STOP 
ENDIF 

IF ( ABS ( DNEW ) . LT . TOLERANCE ) THEN 
SENEW=SESONEW*SO 
ELSE 

SESO=SESONEW 
D=DNEW 
III=III+1 
GOTO  300 
ENDIF 
C 

C  Stress  tensor 

C 

IF (EE .NE . ZERO) THEN 
DO  I=1,NDI,1 

STRESS ( I ) =2 . 0/3 . 0*SENEW/EE*EDEV { I ) 

Sc  +EMOD/ (1.0-2. 0*RNU) ‘EV 

ENDDO 

DO  Iel,NSHR,l 

STRESS ( I+NDI ) =2 . 0/3 . 0*SENEW/EE*EDEV ( I+NDI ) /2 . 0 
ENDDO 
ENDIF 
C 

C  Power law  jacobian 

C 

F=EMOD/ ( 1 . O+RNU+3 . 0/2 . 0* (SENEW/SO) ** (RN-1 . 0 ) ) 

G=EMOD/3 . 0* ( 1 . 0 /( 1 . 0-2 . 0*RNU) -1 . 0 /( 1 . 0+RNU 
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£=  +3.0/2.0*(SENEW/SO)**{RN-1.0)  )  ) 

IF (EE.NE. ZERO) THEN 

H=2 .0/3 . 0*EMOD/EE/ (1 . O+RNU+3 .0/2.0* (SENEW/SO) ** (RN-l . 0) ) 
&  * (1.0/(EE+(RN+1.0)* (SENEW/SO)** (RN-l. 0) 

&  *(SENEW/EMOD) ) -1.0/EE) 

ELSE 

H=ZERO 

ENDIF 

C 

C  Calculate  first  part  of  jacobian  while  also  changing 

C  from  engineering  to  tensor  strains 

C 

DO  I=1,NDI+NSHR,1 
DO  J=l,NDI+NSHR.l 

IF( (I.LE.NDI) .AND. ( J.LE.NDI) )THEN 
DDSDDE  ( I .  J )  =:H*EDEV  ( I )  *EDEV  ( J ) 

ELSEIF( (I.LE.NDI) .AND. (J.GT.NDI) )THEN 
DDSDDE ( I. J)=H*EDEV(I)*EDEV(J) /2.0 
ELSEIF( (I.GT.NDI) .AND. (J.LE.NDI) )THEN 
DDSDDEd,  J)=H*EDEV(I)*EDEV(J)  /2.0 
ELSEIF( (I.GT.NDI) .AND. (J.GT.NDI) )THEN 
DDSDDEd,  J)=H*EDEV(I)*EDEV(J)/2. 0/2.0 
ELSE 

WRITE(7, *) 'Did  not  add  first  jacobian' 

ENDIF 

ENDDO 

ENDDO 

C 

C  Calculate  second  part  of  jacobian 

C 

DO  I=1,NDI,1 
DO  J=1,NDI,1 

DDSDDE (I , J ) =DDSDDE (I , J ) +G 
ENDDO 
ENDDO 
C 

C  Calculate  third  part  of  jacobian 

C 

DO  I=1,NDI+NSHR, 1 

DDSDDE ( I . I ) =DDSDDE (I ,  I ) +F 
ENDDO 
C 

C  Change  jacobian  to  engineering  strain  (2eps  =  gamma) 

C 

DO  I=1,NDI+NSHR,1 
DO  J=1,NDI+NSHR,1 

IF( (I.GT.NDI) .OR. (J.GT.NDI) )THEN 
DDSDDE (I , J ) =DDSDDE (I , J ) / 2 . 0 
ENDIF 
ENDDO 
ENDDO 
C 

C  Update  state  variables 

C 

STATEV ( 1 ) =EE 


STATEV(2)=SENEW 
STATEV ( 3 ) =EE*SENEW 


C 

C  Update  strain  energy,  plastic  dissipation,  and  creep 

C  dissipation 

C 

PRESS=ZERO 
DO  I=1,NDI,1 

PRESS=PRESS+STRESS { I ) /3 . 0 
ENDDO 

SS_1=RK1*S0 

SS_2=RK2*SO 

THETA1=ASIN{ (SS_1/S0-QNC) /R2) 

THETA2=ASIN( (SS_2/S0-QNC) /R2) 

SSE_PRESS=3 . 0* ( 1 . 0-2 . 0*RNU) /2 . 0*PRESS**2 . 0 
SSE_ELASTIC= ( 1 . 0+RNU) /3 . 0*SS_1**2 . 0 
SSE_TRANS= {Rl*R2*SO**2 . 0* (1 . 0/2 . 0*THETA2-1 . 0/4 . 0 
Sc  *SIN(2.0*THETA2) 

Sc  -QNC/R2*  (1.0-{  {SS_2/SO-QNC)  /R2 )  *»2 . 0 )  **0 . 5 )  ) 

&  -(R1*R2*SO**2.0*(1. 0/2. 0*THETA1-1. 0/4.0 

Sc  *SIN{2.0*THETA1) 

Sc  -QNC/R2*{1.0-(  (SS_1/SO-QNC)/R2)**2.0)**0.5)  ) 

SSE_POWER= (( 1 . 0+RNU) /3 . 0*SENEW**2 . 0 
Sc  +RN*SO**2.0/{RN+1.0)*(SENEW/SO)**(RN+1.0)  ) 

Sc  -( (1. 0+RNU) /3.0*SS_2**2.0 

fit  +RN*SO**2.0/(RN+1,0)*{SS_2/SO)**(RN+1.0) ) 

SSE= (SSE_PRESS+SSE_ELASTIC+SSE_TRANS+SSE_POWER) /EMOD 
SPD=ZERO 
SCD=ZERO 
ENDIF 

RETURN 


APPENDIX  B 

CONSTANTS  FOR  SIXTH  ORDER 
RUNGE-KUTTA-VERNER  INTEGRATION 
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The  functions  fji )  are  the  derivatives  of  Yj(6).  The  value  h  in  these 
constants  is  the  step  size  for  integration. 
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APPENDIX  C 

FORTRAN  PROGRAM  FOR  SOLUTION  OF  NONLINEAR 
DIFFERENTIAL  EQUATIONS  WITH  A  SHOOTING  TECHNIQUE 


S’  If’  If  S’ 
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Q  ★★***•*★♦*★♦★★★★★★★**★★**★★**★★★★****★*★★★★♦**♦★♦♦★*****★*★★**★ 

PROGRAM  SHOOTTWO_FIRSTORDER 

^  ************************************************************** 

C 

IMPLICIT  NONE 

CHARACTER*3  CERR, PRI 
INTEGER  I , J , K , ITERATION 

REAL*8  ONE, PIE, PIE2, ZERO, STEPl, STEP2, TOL, ERRF 

REAL* 8  DELTA, PERCENT 

REAL*8  GUESS(5) ,yL (5) , YR (5) ,  PARAM( 50 ) , FEND ( 5) , ENDL ( 5) 

REAL* 8  FENDINIT{5) ,DELTAGUESS(5) ,ENDR(5) 

REAL*8  JAC(5,5),INVJAC{5,5) 

C 

INTEGER  NHARD, COUNT, lUl , IU2 ,  IU3 , IU4 , IU5 , IU6 , IU7 

REAL*8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIAL1(7,2000) ,CONSTANTl{3, 4000) 

COMMON/MAINC/Sl ,  RLAMBDA,  ALPHA,  BETA,  NHARD,  COUNT,  MAXMIS 
COMMON/MAINA/POTENTIALl , CONSTANTl 
COMMON/OUT/ lUl, IU2, IU3, 1U4, IU5, IU6, IU7 
C 

CALL  INPUT  (TOL, PARAM, STEPl , STEP2 , ERRF , DELTA, NHARD, 

&  RLAMBDA, PERCENT, GUESS ) 


C 


SI  =  -1 . 0/ (NHARD+1)  !  First  eigen  value 

ALPHA  =  RLAMBDA** 2 -RLAMBDA+1 

BETA  =  2*RLAMBDA**2-2*RLAMBDA-1 

ZERO  =  0.0 

ONE  =  1.0 

PIE  =  4*DATAN(ONE) 

PIE2  *  PERCENT* PIE 

COUNT  =  0 

MAXMIS  =  -1000000.0 


ITERATION  =  1 
CALL  OPENUNIT 
10  CALL  FIRSTTRY 


CALL  NORMALIZE 

CALL  SORT 

CALL  OOTPUT 

CALL  SORT 

CALL  INTEGCONST 

IF(CERR(1:1) .EQ. ’N* ) 


{ ITERATION, lUl , IU2 , IU3 , IU4 , IU5 , IU6 , IU7 ) 
(CERR, PRI , ITERATION, PIE, PIE2 , ZERO, STEPl , 
STEP2 , TOL , ERRF , GUESS , YL , YR , PARAM , FEND , 
ENDL. FENDINIT, ENDR) 

(COUNT, MAXMIS , POTENTIALl ) 

( COUNT, POTENTl ALl , 7 , 2 0 0 0 ) 

( 2 *COUNT, CONSTANTl , 3 , 4 0 0 0 ) 

( 2 *COUNT, CONSTANTl ) 

GOTO  20 


IMPROVE  (CERR, PRI, ITERATION, PIE, PIE2, ZERO, STEPl, 

STEP2 , TOL , ERRF , CUES S , YL , YR , PARAM , FEND , 
ENDL, FENDINIT, ENDR, DELTA, JAC) 

CALL  INVERTJACOBIAN  (JAC.INVJAC) 

CALL  FINDDELTAGUESS  ( INVJAC , FENDINIT , DELTAGUESS ) 

CALL  ADDGUESSANDDELTA  (GUESS, DELTAGUESS) 

CALL  CLOSEUNIT  (lUl, IU2, IU3, IU4, IU5, IU6, IU7) 

CALL  OPENUNIT  (1, lUl, IU2, IU3 , IU4, IU5, IU6, IU7) 

ITERATION  =  ITERATION+1 
COUNT=0 
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c 

c 


c 

c 


c 

c 


c 

c 


c 

c 


c 

c 


c 

c 


MAXMIS=-1000000.0 
GOTO  10 

20  CALL  CLOSEUNIT  (lUl, IU2, IU3, IU4, 1U5, IU6, 1U7) 

STOP 

END 

ii************************************************************* 

SUBROUTINE  INPUT (TOL, PARAM, STEPl, STEP2 , ERRF, DELTA, NHARD, 

&  RLAMBDA , PERCENT , GUESS ) 

**********************'********'**«***»********'**********«****** 


IMPLICIT 

INTEGER 

REAL*8 

REAL*8 


NONE 

NHARD 

STEPl , STEP2 , TOL, ERRF, DELTA, PERCENT, RLAMBDA 
PARAM{50) ,GUESS{5) 


OPEN  (UNIT=1 , FILE= • INPUT . DAT ‘ , STATUS= • OLD ‘ ) 


READd,  *) 
READd,  *) 
READ (1 , * ) 
READd,  *) 
READd,  *) 
READd,  *) 
READd,*) 
READd,*) 
READd,*) 
READd,  *) 
READd,  *) 
READd,  *) 
READd,*) 
READd,*) 
READd,  *) 
READd,*) 
READd,*) 


TOL 

PARAMd) 

PARAM(2) 

PARAMO) 

PARAMdO) 

STEPl 

STEP2 

ERRF 

DELTA 

NHARD 

RLAMBDA 

PERCENT 

GUESS ( 1 ) 

GUESS (2) 

GUESS (3) 

GUESS (4) 

GUESS (5) 


Global  error  tolerence  <  0,001 

Initial  guess  of  step  size 

Absolute  value  of  minimum  step  size 

Absolute  value  of  maximum  step  size 

Type  of  error  control 

Number  of  increments  for  zero  to  pie/2 

Number  of  increments  for  pie  to  pie/2 

Acceptable  difference  at  pie/2 

Increment  for  guesses 

Ramberg-Osgood  hardening  coefficient 

Measure  of  the  degree  of  plasticity 

Percent  of  pie  to  shoot  too 

Five  initial  guess  for  boundary  conditions 


CLOSE  (UNIT=1) 


RETURN 

END 


SUBROUTINE  FIRSTTRY  (CERR, PRI, ITERATION, PIE, PIE2 , ZERO, STEPl , 

&  STEP2, TOL, ERRF, GUESS, YL, YR, PARAM, FEND, 

Sc  ENDL ,  FENDINIT ,  ENDR ) 

♦************************************************************* 


IMPLICIT  NONE 

CHARACTER*3  CERR, PRI 
INTEGER  ITERATION 

REAL* 8  PIE, PIE2, ZERO, STEPl, STEP2, TOL, ERRF 

REAL*8  GUESS(5) , YL(5) ,YR(5) , PARAM(50) ,FEND(5) ,ENDL(5) 

REAL*8  FENDINIT(5) ,ENDR(5) 
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c 

PRI  =  'YES' 

CALL  LEFTBOUND  ( GUESS , YL ) 

CALL  RIGHTBOUND  (GUESS, YR) 

CALL  im’EGRATZEROPlE2  (ZERO, P1E2 , STEPl , TOL, PARAM, YL, PRI ) 

CALL  INTEGRATPIEPIE2  (PIE, PIE2, STEP2, TOL, PARAM, YR, PRI) 

CALL  SAVEENDPOINTS  (YL, YR, ENDL, ENDR) 

CALL  FENDPOINTS  (YL, YR, FENDINIT) 

CALL  FERROR  (ITERATION, FENDINIT, ERRF, CERR) 

C 

RETURN 

END 

C 

Q  ****ii-k*tr**1i********ii**-k**it***ir*ift*it**1r**ir*it*****1i*****ii****ii** 

SUBROUTINE  IMPROVE  (CERR, PRI, ITERATION, PIE, PIE2 , ZERO, STEPl, 

Sc  STEP2, TOL, ERRF, GUESS,  YL,YR,  PARAM, FEND, 

&  ENDL , FENDINIT , ENDR , DELTA , JAC ) 

Q  it************************************************************* 

c 

IMPLICIT  NONE 

CHARACTER*!  CERR, PRI 
INTEGER  I, ITERATION 

REAL*8  PIE, PIE2 , ZERO, STEPl , STEP2 , TOL, ERRF, DELTA 

REAL*8  GUESS (5) ,YL(5) , YR(5) , PARAM(50) ,FEND(5) ,ENDL(5) 

REAL*8  FENDINIT (5) ,ENDR(5) 

REAL*8  JAC(5,5) 

C 

PRI  =  'NOO' 

DO  1=1, 5,1 

CALL  INCREMENTGUESS  ( DELTA, GUESS ( I ) ) 

CALL  LEFTBOUND  (GUESS, YL) 

CALL  RIGHTBOUND  (GUESS, YR) 

CALL  INTEGRATZEROPIE2  (ZERO, PIE2 , STEPl , TOL, PARAM, YL, PRI) 
CALL  INTEGRATPIEPIE2  (PIE, P1E2 , STEP2 , TOL, PARAM, YR, PRI) 

CALL  FENDPOINTS  (YL,yR,FEND) 

CALL  JACOBIAN  (I, DELTA, FENDINIT, FEND, JAC) 

CALL  UNINCREMENTGUESS  (DELTA. GUESS ( I ) ) 

ENDDO 

C 

RETURN 

END 

C 

Q  ********************ir*1rii****1r**t,**1r************ir*************it 

SUBROUTINE  LEFTBOUND  (GUESS, YL) 

C  ************************************************************** 

c 

IMPLICIT  NONE 

INTEGER  I 

REAL*8  GUESS(5) ,YL(5) 

C 

WRITEdl,*)  'Setting  left  boundary  conditions' 

C 

YL(1)  =  1.0 

YL(2)  =  0 
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YL(3)  =  GUESS(l) 

YL(4)  =  0.0 

YL(5)  =  GUESS(2) 

C 

RETURN 

END 

C 

Q  ******'******'***********•*'***************************«********* 

SUBROUTINE  RIGHTBOUND  (GUESS, YR) 

Q  *************************************************************** 

C 

IMPLICIT  NONE 

INTEGER  I 

REALMS  GUESS(5) ,YR(5) 

C 

WRITE (11,*)  'Setting  right  boundary  conditions' 

C 

YR ( 1 )  =  0 

YR(2)  =  0 

YR(3)  =  GUESS (3) 

YR(4)  =  GUESS (4) 

YR(5)  =  GUESS(5) 

C 

RETURN 

END 

C 

C  ************************************************************** 

SUBROUTINE  INTEGRATZEROPIE2  (ZERO, PIE2, STEPl, TOL, PARAM, YL, PRI) 

C  «*«*«'***«*******«********««*****«**«***«**«******«*********'**i^ 

c 

IMPLICIT  NONE 
CHARACTER*3  PRI 
INTEGER  I 

REAL* 8  ZERO, PIE2 , STEPl , X, TOL, RINCREMENTl , OVER 

REAL*8  YL(5) ,PARAM(50) ,DYDX(5) 

C 

WRITE (11,*)  'Integrating  from  zero  to  pie/2' 

C 

RINCREMENTl  =  (PIE2-ZERO) /INT(STEPl) 

IF  (PRI(1:1) .EQ. 'Y' )  THEN 
CALL  SAVE  (ZERO,YL) 

ENDIF 

OVER  =  RINCREMENTl/ 2 

DO  XsZERO, PIE2-RINCREMENT1+OVER, RINCREMENTl 

CALL  INTEGRATE  (X, X+RINCREMENTl, RINCREMENTl , TOL, PARAM, YL) 

IF  (PRI(1;1) .EQ. 'Y' )  THEN 

CALL  SAVE(X+RINCREMENT1,YL) 

ENDIF 

ENDDO 

C 

RETURN 

END 

C 
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c 

c 

c 


c 

c 


c 


c 

c 

c 

c 


c 


c 


c 


c 

c 

c 

c 


************************************************************** 

SUBROUTINE  INTEGRATPIEPIE2  ( PIE, PIE2 , STEP2 , TOL, PARAM, YR, PRI) 

************************************************************** 


IMPLICIT  NONE 
CHARACTER*3  PRI 
INTEGER  I 

REAL*8  PIE, PIE2 , STEP2, X, TOL, RINCREMENT2 , OVER 

REAL*8  YR{5) ,PARAM(50) ,DyDX(5) 

WRITE (11,*)  ‘Integrating  from  pie  to  pie/2' 

RINCREMENT2  =  (PIE2-PIE) /INT(STEP2) 

IF  (PRI(1:1) .EQ. ‘Y* )  THEN 
CALL  SAVE ( PIE, YR) 

ENDIF 

OVER  =  RINCREMENT2/2 

DO  X=PIE, PIE2-RINCREMENT2+OVER,RINCREMENT2 

CALL  INTEGRATE  (X, X+RINCREMENT2, RINCREMENT2 , TOL, PARAM, YR) 

IF  (PRI{1:1) .EQ. ‘Y' )  THEN 

CALL  SAVE{X+RINCREMENT2,YR) 

ENDIF 

ENDDO 

RETURN 

END 

************************************************************** 

SUBROUTINE  INTEGRATE  (BEGIN, END, STEP, TOL, PARAM, Y) 

************************************************************** 


IMPLICIT 

INTEGER 

REAL*8 

REAL*8 

EXTERNAL 


NONE 

IDO,NEQ 

BEGIN, END, STEP, XEND,  X, TOL,  OVER 
Y  ( 5 ) ,  Y  JUNK  ( 5 ) ,  FCN ,  P ARAM  ( 5  0 ) 
FCN,DIVPRK 


IDO  =  1 
NEQ  =  5 
OVER  =  STEP/ 2 

DO  X=BEGIN, END-STEP+OVER, STEP 
XEND  =  X+STEP 

CALL  DIVPRK  ( IDO, NEQ, FCN, X, XEND, TOL, PARAM, Y) 
ENDDO 


IDO  =  3 

CALL  DIVPRK  ( IDO, NEQ, FCN, X, XEND, TOL, PARAM, Y JUNK) 

RETURN 

END 


************************************************************** 

SUBROUTINE  SAVEENDPOINTS  (YL, YR, ENDL, ENDR) 

************************************************************** 
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c 

c 


c 


c 

c 

c 

c 


c 

c 


c 


c 

c 

c 

c 


c 


c 

c 


IMPLICIT  NONE 

INTEGER  I 

REAL*8  YL(5) ,YR(5) ,ENDL(5) ,ENDR(5) 

WRITE (11,*)  'Saving  endpoints' 

DO  1=1, 5,1 

ENDL(I)  =  YL(I) 

ENDR(I)  =  YR{I) 

ENDDO 

RETURN 

END 

♦************************************************************* 

SUBROUTINE  FENDPOINTS  (YL,YR,FEND) 

************************************************************** 


IMPLICIT  NONE 
INTEGER  I 

REAL*8  YL{5) ,YR(5) ,FEND{5) 

WRITE (11,*)  'Determining  error  in  the  endpoints' 
DO  1=1, 5,1 

FEND(I)  =  YR(I)-YL(I) 

ENDDO 

RETURN 

END 


♦♦It*********************************************************** 

SUBROUTINE  FERROR  ( ITERATION, FEND, ERRF, CERR) 

************************************************************** 


IMPLICIT  NONE 

CHARACTER*3  CERR 
INTEGER  I, ITERATION 
REAL*8  ERRMAX, ERRF, ERR 

REAL*8  FEND(5) 

INTEGER  lUl , IU2 , IU3 , IU4 , IU5 , IU6 , IU7 
COMMON/ OUT/ lUl, IU2, IU3 , IU4, IU5, IU6, IU7 

WRITE (11,*)  'Determining  maximimi  error  in  F  at  pie/2' 

ERRMAX  =  0 
DO  1=1, 5,1 

ERR  =  DABS(FEND(I) ) 

ERRMAX  =  MAX (ERR, ERRMAX) 

ENDDO 

WRITE(*,*)  'Max  error  for  iteration  ', ITERATION, '  is  ', ERRMAX 
WRITEdl,*)  'Max  error  for  iteration  ',  ITERATION,  '  is  ', ERRMAX 
IF  ( ERRMAX. LE. ERRF)  THEN 
CERR  =  'NOO' 


165 


c 


c 

c 

c 

c 


c 

c 


c 

c 

c 

c 


c 

c 


c 

c 

c 

c 

c 


c 


.c 


c 

c 

c 


ELSE 

CERR  =  'YES' 

ENDIF 

RETURN 

END 

******’******************************************************** 

SUBROUTINE  INCREMENTGUESS  (DELTA, G) 

************************************************************** 

IMPLICIT  NONE 

REAL* 8  DELTA, G 

G  =  G  +  DELTA 

RETURN 

END 

SUBROUTINE  UNINCREMENTGUESS  (DELTA, G) 

************************************************************** 

IMPLICIT  NONE 

REAL* 8  DELTA, G 

G  =  G  -  DELTA 

RETURN 

END 


SUBROUTINE  JACOBIAN  (ICOL, DELTA, FENDINIT, FEND, JAC) 

************************************************************** 


IMPLICIT 

INTEGER 

REAL*8 

REAL*8 

REAL*8 


NONE 
I, ICOL 
DELTA 

FENDINIT ( 5 ), FEND ( 5 ) 
JAC(5,5) 


DO  1=1, 5,1 

JAC ( I , ICOL ) = ( FEND ( I ) -FENDINIT ( I ) ) /DELTA 
ENDDO 


RETURN 

END 


SUBROUTINE  INVERT JACOB IAN  (JAC,INVJAC) 

•k************************************************************* 


c 

c 


IMPLICIT 


NONE 
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c 

c 


c 

c 

c 

c 


c 

c 


c 


c 

c 

c 

c 


c 

c 


c 


c 

c 

c 

c 


INTEGER  I ,  J 

REAL*8  JAC(5.5) ,INVJAC(5,5) 

CALL  DLINRG  (5, JAC, 5, INVJAC, 5) 

RETURN 

END 

SUBROUTINE  FINDDELTAGUESS  ( INVJAC, FENDINIT, DELTAGUESS ) 

************************************************************** 


IMPLICIT 

INTEGER 

REAL*8 

REAL*8 


NONE 
I,  J 

FENDINIT { 5 ) , DELTAGUESS { 5 ) 
INVJAC (5, 5) 


V7RITE(11,*)  ‘Determine  change  in  guesses' 

DO  1=1, 5,1 

DELTAGUESS { I )  =  0 
DO  J=l,5,l 

DELTAGUESS (I)  =  DELTAGUESS ( I) -INVJAC ( I, J) ‘FENDINIT (J) 
ENDDO 
ENDDO 


RETURN 

END 

*************************************************************** 

SUBROUTINE  ADDGUESSANDDELTA  (GUESS, DELTAGUESS) 

***•*«***«*•**«******•*****•*•***«****««***•**•***'*****•»***** 


IMPLICIT  NONE 

INTEGER  I 

REAL*  8  GUESS ( 5 ) , DELTAGUESS ( 5 ) 


WRITE (11,*)  ‘Adding  the  change  in  guesses  to  previous  guesses' 


GUESS { 1 ) 
GUESS (2) 
GUESS (3) 
GUESS (4) 
GUESS (5) 


GUESS ( 1 ) +DELTAGUESS ( 1 ) 
GUESS ( 2 ) +DELTAGUESS ( 2 ) 
GUESS ( 3 ) +DELTAGUESS ( 3 ) 
GUESS ( 4 ) +DELTAGUESS ( 4 ) 
GUESS ( 5 ) +DELTAGUESS ( 5 ) 


RETURN 

END 


SUBROUTINE  FCN{NEQ,XX,Y,DYDX) 


IMPLICIT  NONE 

INTEGER  NEQ,  I 

REAL*8  XX,Y{NEQ) ,DYDX{NEQ) ,MEAN 
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REAL*8  SIG(4) , DSIG (4 ) , DDSIG ( 4 ) , S (4 ) ,MIS(5) 

REAL*8  C11,C12 

REAL* 8  CN11,CN12,CN13,CN14,CN15 

REAL*8  CL11,CL12,CL13,CL14,CL15 

REAL*8  C21,C22,C23,C24.C25 

C 

INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIAL1(7,2000) ,CONSTANT1(3,4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/ POTENTIALl , CONSTANTl 
C 

C  Find  the  values  of  the  last  field  at  X 
C 

CALL  STRESS(Y(5) ,Y,SIG) 

CALL  DSTRESS{Y{5) ,Y,DSIG) 

CALL  DDSTRESS(Y(5) ,Y,DDS1G) 

CALL  DEV(SIG,S,MEAN) 

CALL  SIGMIS (SIG, SIG, DSIG, DSIG, DDSIG, DDSIG, MIS) 
MIS(1)=MIS(1)**0.5 
C 

Cll  =  2*ALPHA*(DSIG(1)**2+SIG(2)*DDSIG(2)+DSIG(2)**2) 
&  +  BETA*{SIG(1)*DDSIG{2)+2*DSIG(1)*DSIG(2) ) 

£c  +  6*(SIG{3)*DDSIG{3)+DSIG{3)**2) 

t  +  (Y(5)+2)*{2*ALPHA*SIG(1)+BETA*SIG(2) )*Y(3) 

C12  =  2*ALPHA*SIG(1)+BETA*SIG{2) 

C 

C  Linear  coefficients  of  the  differential  equation 
C 

CLll  =  NHARD*Y(5)*{((l+RLAMBDA)/3+(NHARD*Y(5)+l) 

Sc  *(2-RLAMBDA)/3)*(Y(5)+l)-(  (2-RLAMBDA)/3 

&  +(NHARD*Y(5)+l)*(l+RLAMBDA)/3) )*{y(5)+2) 

CL12  =0.0 

CL13  =  (2-RLAMBDA)/3*(Y{5)+2) 

&  -(1+RLAMBDA)/3*{Y(5)+2)*(Y(5)+1) 

&  +2*(NHARD*Y(5)+1)*(Y(5)+1) 

&  -NHARD*Y(5)*{ {2-RLAMBDA)/3+(NHARD*Y{5)+l) 

&  *(l+RLAMBDA)/3) 

CL14  =0.0 
CL15  =  (2-RLAMBDA)/3 
C 

C  Nonlinear  coefficients  of  the  differential  equation 
C 

CNll  =  (NHARD-l)*(Y(5)+2)*MIS(l)**(-2) 

&  *(Cll/2+(NHARD-3)*MIS(3)**2) 

&  *( (2-RLAMBDA)/3-(l+RLAMBDA)/3*(Y(5)+l) ) 

CN12  =  2*(NHARD-1)*MIS(1)**(-1)*MIS(3) 

&  *{ (2-RLAMBDA)/3*(Y(5)+2) 

&  -(1+RLAMBDA)/3*(Y{5)+2)*(Y(5)+1) 

Sc  +(NHARD*Y(5)+1)*(Y(5)+1)  ) 

CN13  =  (2-RLAMBDA)/3*(NHARD-l)*MIS(l)**(-2) 

Sc  *(Cll/2+{NHARD-3)*MIS(3)**2) 

CN14  =  2* (NHARD-l) * {2-RLAMBDA) /3*MIS(1) ** (-1) *MIS{3) 
CN15  =  C12*(NHARD-1)/2*MIS(1)**(-2)*S{1) 


C 
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C  Coefficients  for  the  differential  equation 
C 

C21  =  CLll+CNll 
C22  =  CL12+CN12 
C23  =  CL13+CN13 
C24  =  CL14+CN14 
C25  =  CL15+CN15 
C 

C  Set  of  first  order  differential  equations 

C 

DYDXd)  =  y(2) 

DYDX(2)  =  Y(3) 

DYDX(3)  =  Y{4) 

DYDX{4)  =  -(C24*Y(4)+C23*Y(3)+C22*Y(2)+C21*Y(1) ) /C25 
DYDX(5)  =0.0 
C 

RETURN 

END 

C 

SUBROUTINE  STRESS (SI , F, SIG) 

C  ************************************************************** 

c 

IMPLICIT  NONE 

REAL*8  SI 

REAL*8  F{5),SIG(4) 

C 

INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA, BETA,  MAXMIS 

REAL*8  POTENTIALl (7 , 2000) , CONSTANTl (3,4000) 

COMMON/MAINC/Sl ,  RLAMBDA,  ALPHA,  BETA,  NHARD,  COUNT,  MAXMIS 
COMMON/MAINA/ POTENTIALl , CONSTANTl 
C 

C  Determine  the  stresses 
C 

SIG(1)=F(3)+(SI+2)*F(1) 

SIG(2)=(SI+2)*(SI+1)*F(1) 

SIG(3)=-(SI+1)*F(2) 

SIG(4)=RLAMBDA*(SIG(1)+SIG(2) ) 

C 

RETURN 

END 

C 

SUBROUTINE  DSTRESS (SI, F, DSIG) 

C  ************************************************************** 

c 

IMPLICIT  NONE 

REAL*8  SI 

REAL*8  F(5),DSIG(4) 

C 

INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIALl (7, 2000) .CONSTANTl (3,4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
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COMMON/MAINA/ POTENTIALl , CONSTANTl 
C 

C  Determine  the  first  derivative  of  the  stresses 
C 

DS1G(1)=F{4)+{SI+2)*F(2) 

DSIG{2)=(SI+2)*{SI+1)*F(2) 

DSIG(a)=-(SI+l)*F{3) 

DSIG(4)=RLAMBDA*(DSIG(1)+DSIG(2) ) 

C 

RETURN 

END 

C 

Q  ************1r***************1fit**itit************¥t************‘k** 

SUBROUTINE  DDSTRESS ( S I , F , DDS IG ) 

Q  ***1tit**1r*******iHt***********iflt******it*1t*****irir1t**ir**'lt**ir***1r** 

C 

IMPLICIT  NONE 

REAL*8  SI 

REAL*8  F(5) ,DDSIG(4) 

C 

INTEGER  NHARD, COUNT 

REAL»8  RLAMBDA, SI. ALPHA, BETA, MAXMIS 

REAL*8  POTENTIALl (7, 2000) .CONSTANTl (3,4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/ POTENTIALl , CONSTANTl 
C 

C  Determine  the  second  derivative  of  the  stresses 

C 

DDSIG(1)=F{5)+(SI+2)*F(3) 

DDSIG(2)=(SI+2)*{SI+1)*F(3) 

DDSIG{3)=-(SI+1)*F(4) 

DDSIG { 4 ) =RLAMBDA* (DDSIG ( 1 ) +DDSIG ( 2 ) ) 

C 

RETURN 

END 

C 

^  ****ir******irir*ir*ir***»**ir**1rit***1rint**********1c*****r****1r****1rit1r 

SUBROUTINE  DEV ( S IG, S, MEAN) 

Q  ***************ift****************************ii*1Hr**ir****1i***** 

C 

IMPLICIT  NONE 

REAL*8  MEAN 

REAL*8  SIG(4),S(4) 

C 

INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIALl (7,2000) , CONSTANTl (3,4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/ POTENTIALl , CONSTANTl 
C 

C  Determine  the  deviatoric  stresses 

C 

S(l)  =  (2-RLAMBDA)/3*SIG(l)-(l+RLAMBDA)/3*SIG(2) 

S(2)  =  (2-RLAMBDA)/3*SIG(2)-(l+RLAMBDA)/3*SIG(l) 

S(3)  =  S1G(3) 
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c 


c 

c 

c 

c 


c 


c 

c 

c 


c 


c 

c 

c 

c 


c 


c 

c 

c 


c 


S{4)  =  (2*RLAMBDA-1)/3*(SIG(1)+SIG{2)) 
MEAN  =  (1+RLAMBDA)/3*(SIG(1)+SIG(2)) 

RETURN 

END 


SUBROUTINE  DDEV(DSIG, DS) 


IMPLICIT  NONE 

REAL*8  DSIG(4) .DS(4) 


INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIAL1(7,2000) ,CONSTANTl (3,4000) 

COMMON/MAINC/Sl,  RLAMBDA,  ALPHA,  BETA,  NHARD,  COUNT, MAXMIS 
COMMON/MAINA/ POTENTIALl , CONSTANTl 


Determine  the  first  derivative  of  the  deviatoric  stresses 

DS(1)  =  (2-RLAMBDA)/3*DSIG(l)-(l+RLAMBDA)/3*DSIG{2) 

DS(2)  =  (2-RLAMBDA) /3*DSIG{2)-(1+RLAMBDA) /3*DSIG(1) 

DS(3)  =  DSIG(3) 

DS(4)  =  (2*RLAMBDA-1)/3*(DSIG(1)+DSIG(2) ) 

RETURN 

END 

************************************************************** 

SUBROUTINE  DDDEV { DDS IG , DDS ) 

*************************************************************** 


IMPLICIT  NONE 

REAL*8  DDSIG{4) ,DDS(4) 


INTEGER  NHARD, COUNT 

REAL* 8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIALl (7, 2000 ), CONSTANTl (3, 4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/ POTENTIALl , CONSTANTl 


Determine  the  second  derivative  of  the  deviatoric  stresses 


DDS(l)  =  (2-RLAMBDA) /3*DDSIG(1)-(1+RLAMBDA)/3*DDSIG(2) 
DDS(2)  =  (2-RLAMBDA) /3*DDSIG(2)-(1+RLAMBDA)/3*DDSIG(1) 
DDS (3)  =  DDSIG(3) 

DDS (4)  =  (2*RLAMBDA-1) /3* (DDSIG{1)+DDSIG(2) ) 


C 


RETURN 

END 
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c 

c 

c 


c 


c 

c 

c 


c 


c 

c 

c 

c 


c 


************************************************************** 

SUBROUTINE  SIGMIS (SIGl , SIG2 , DSIGl , DSIG2 , DDSIGl , DDSIG2 , MIS ) 

************************************************************** 


IMPLICIT  NONE 

REAL*8  SIG1(4),DSIG1(4),DDSIG1{4) ,MIS(5) 

REAL*8  SIG2(4) ,DSIG2(4) ,DDSIG2(4) 


INTEGER  NHARD, COUNT 

REAL* 8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIALl{7, 2000) ,CONSTANTl(3, 4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/POTENTIALl , CONSTANTl 


Determine  the  angular  functions  for  the  effective  stress 


MIS(l)  =  ALPHA*(SIG1(1)*SIG2{1)+SIG1(2)*SIG2(2) ) 

&  +BETA/2*(SIG1{1)*SIG2(2)+SIG2(1)*SIG1{2) ) 

&  +3*SIG1(3)*SIG2(3) 

MIS(2)  =  ALPHA*(SIG1(1)*DSIG2(1)+DSIG1{1)*SIG2(1) 

&  +SIG1(2)*DSIG2(2)+DSIG1(2)*SIG2(2) ) 

St  +BETA/2*  (SIGl(l)  *DSIG2  (2)+DSIGl  (1)  *SIG2  (2) 

St  +SIG2{1)*DSIG1(2)+DSIG2(1)*SIG1{2)  ) 

St  +3*  (SIG1(3)*DSIG2(3)+DSIG1(3)*SIG2(3)  ) 

1) ,GT.0)THEN 
)  =  MIS(2)/2/MIS{l)**0.5 


IF (MIS ( 
MIS  (3 
ENDIF 
MIS (4) 


St 

St 

St 

St 

St 

St 

Sc 

Sc 

Sc 


IF (MIS ( 
MIS  (5 
ENDIF 


=  ALPHA* ( SIGl ( 1 ) *DDS1G2 ( 1 ) +DSIG1 ( 1 ) *DSIG2 ( 1 ) 
+DSIG1 ( 1 ) *DSIG2 ( 1 ) +DDSIG1 ( 1 ) *SIG2 ( 1 ) 
+SIG1(2)*DDSIG2(2)+DSIG1(2)*DSIG2(2) 
+DS IGl { 2 ) *DSIG2 { 2 ) +DDS IGl ( 2 ) *  S IG2 { 2 ) ) 
+BETA/2* (SIGl (1) *DDSIG2 (2 ) +DSIG1 ( 1 ) *DSIG2 (2) 
+DSIG1 ( 1 ) *DSIG2 ( 2 ) +DDSIG1 ( 1 ) *SIG2 ( 2 ) 
+SIG2 ( 1 ) *DDSIG1 ( 2 ) +DSIG2 ( 1 ) *DSIG1 ( 2 ) 
+DSIG2(1)*DSIG1(2)+DDSIG2(1)*SIG1(2) ) 
+3* (SIGl (3) *DDSIG2 (3)+DSIGl (3 ) *DSIG2 (3 ) 
+DSIG1(3)*DSIG2(3)+DDSIG1(3) *SIG2(3) ) 

1) .GT.0)THEN 

)  =  (MIS(4)-MIS(2)**2/2/MIS(l) )/2/MIS(l)**0.5 


RETURN 

END 

************************************************************** 
SUBROUTINE  STRAIN (MIS, S, EPS) 

************************************************************** 


IMPLICIT  NONE 

REAL*8  MIS(5) ,S(4) ,EPS(4) 


INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL  *  8  POTENTI ALl (7,2000), CONSTANTl (3,4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
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COMMON/MAINA/POTENTIALl , CONSTANTl 
C 

C  Determine  the  angular  strain  functions 
C 

EPS(l)  =  3*MIS(l)**(NHARD-l)*S(l)/2 
EPS(2)  =  3*MIS(l)**{NHARD-l)*S{2)/2 
EPS(3).  =  3*MIS(l)**(NHARD-l)*S(3)/2 
EPS (4)  =  0 
C 

RETURN 

END 

C 

0  ***********«***********•*'*******•'***************«*********'**'** 

SUBROUTINE  DSTRAIN (MIS, S, DS, DEPS) 

Q  it************************************************************* 

C 

IMPLICIT  NONE 

REAL*8  MIS(5) ,S(4)  ,DS(4)  ,DEPS{4) 

C 

INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIALl (7, 2000) .CONSTANTl (3, 4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA,  BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/POTENTIALl , CONSTANTl 
C 

C  Determine  the  first  derivative  of  the  angular  strain  functions 

C 

DEPS(l)  =  3*(MIS(1)**(NHARD-1)*DS(1)+(NHARD-1)*MIS(1)**{NHARD- 

2) 

&  *MIS(3)*S(1) )/2 

DEPS{2)  =  3*(MIS(1)**(NHARD-1)*DS(2)+(NHARD-1)*MIS{1)**(NHARD- 

2) 

&  *MIS{3)*S(2) )/2 

DEPS(3)  =  3*(MIS(1)**(NHARD-1)*DS(3)+(NHARD-1)*MIS(1)**(NHARD- 

2) 

£c  *MIS(3)*S(3)  )/2 

DEPS(4)  =  0 
C 

RETURN 

END 

C 

^  •Ir**************’**************************************'^******* 

SUBROUTINE  DDSTRAIN (MIS , S , DS , DDS , DDEPS ) 

C  ************************************************************** 

c 

IMPLICIT  NONE 

REAL*8  MIS{5) ,S(4) ,DS(4) ,DDS(4) ,DDEPS(4) 

C 

INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIALl(7, 2000) ,CONSTANTl(3, 4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/POTENTIALl , CONSTANTl 
C 

C  Determine  the  second  derivative  of  the  angular 
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C  strain  functions 

C 

DDEPS(l) 

£c 

Sc 
Sc 
Sc 
Sc 

DDEPS{2) 

Sc 
Sc 
Sc 
Sc 
& 

DDEPSO) 

Sc 
Sc 
Sc 
Sc 
Sc 

DDEPS { 4 ) 

C 

RETURN 
END 
C 

Q  ************************************************************** 

SUBROUTINE  DISPLACE (EIGEN, MIS, EPS, DEPS, U) 

Q  ************************************************************** 

c 

IMPLICIT  NONE 

REAL*8  EIGEN 

REAL*8  U(2) ,EPS(4) ,MIS(5) ,DEPS{4) 

C 

INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*  8  POTENTI ALl (7,2000), CONSTANTl (3,4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/ POTENTIALl , CONSTANTl 
C 

C  Detennine  the  angular  displacement  functions 

C 

U(l)  =  EPS(1)/(NHARD*EIGEN+1) 

U(2)  =  (2*EPS(3)-DEPS(1)/(NHARD*EIGEN+1) )/(NHARD*EIGEN) 

C 

RETURN 

END 

C 

^  ************************************************************** 

SUBROUTINE  DDISPLACE (EIGEN, MIS, EPS, DEPS, DDEPS, DU) 

Q  ************************************************************** 

C 

IMPLICIT  NONE 

REAL*8  EIGEN 

REAL*8  DU(2) ,EPS(4) ,MIS(5) ,DEPS(4) ,DDEPS(4) 

C 


=  3*(MIS(1)**(NHARD-1)*DDS(1)+ 

(NHARDl ) *MIS ( 1 ) ** (NHARD-2 ) 

♦MIS ( 3 ) *DS ( 1 ) + (NHARD-1 ) * (MIS ( 1 ) * * (NHARD-2 ) 
*(MIS(3) *DS(1)+MIS(5) *S(1) )+ 

(NHARD2 ) *MIS ( 1 ) * * (NHARD-3 ) 

*MIS(3)**2*S{1) ))/2 
=  3*(MIS(1)**(NHARD-1)*DDS{2)+ 

(NHARDl) *MIS{1) ♦* (NHARD-2) 

♦MIS ( 3 ) *DS ( 2 )+ (NHARD- 1 )* (MIS ( 1 )**( NHARD-2 ) 
♦(MIS(3) ♦DS(2)+MIS(5)^S(2) )+ 

(NHARD2 ) ♦MIS { 1 ) ♦♦ (NHARD-3 ) 
♦MIS(3)^^2^S{2)))/2 
=  3^(MIS(1)^^(NHARD-1)^DDS(3)+ 

(NHARDl) ♦MIS (1) ♦♦ (NHARD-2) 

♦MIS ( 3 ) ♦DS ( 3 )+ (NHARD- 1 )♦ (MIS ( 1 )♦♦( NHARD-2 ) 
♦(MIS(3)^DS(3)+MIS(5)^S(3) )+ 

(NHARD2 ) ♦MIS ( 1 ) ♦♦ (NHARD-3 ) 
♦MIS(3)^^2^S(3)))/2 
=  0 


174 


c 

c 

c 


c 


c 

c 

c 

c 


c 


c 

c 

c 


c 

c 


c 

c 

c 

c 


INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIAL1(7,2000) ,CONSTANT1(3,4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/POTENTIALl , CONSTANTl 


Determine  the  derivative  of  the  angular  displacement  functions 
DU{1)  =  DEPS(1)/(NHARD*EIGEN+1) 

DU(2)  =  (2*DEPS(3)-DDEPS{1)/(NHARD*E1GEN+1) )/{NHARD*EIGEN) 

RETURN 

END 

*********<***«*«****'»'*«********ll^*'**************t^**t^************ 

SUBROUTINE  INTEGRAND (X, EIGEN, MIS , SIGl , U1 , DUl , KOUNT) 

**************************************************************** 


IMPLICIT 

INTEGER 

REAL*8 

REAL*8 


NONE 
KOUNT 
X, EIGEN 

MIS (5) ,SIG1(4) ,U1(2) ,DU1(2) 


INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIALl (7,2000) , CONSTANTl (3, 4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/POTENTIALl , CONSTANTl 


Integrand  for  the  J- integral 


CONSTANTl ( 1 , KOUNT) =X 

CONSTANTl (2, KOUNT) =NHARD*MIS { 1 ) ** (NHARD+1 ) *DCOS (X) / (NHARD+1) 

&  -DSIN(X)*(SIG1(1)*(U1(2)-DU1(1) )-SIGl{3)*(Ul(l)+DUl(2) ) ) 

Sc  -DCOS  (X)  *  (NHARD*EIGEN+1 )  *  (SIGl  ( 1 )  *U1  ( 1 )  +SIG1  ( 3 )  *U1  (2 )  ) 

CONSTANTl ( 1 , KOUNT+1 ) =-X 

CONSTANTl ( 2 , KOUNT+ 1 ) =NHARD*MIS (!)**( NHARD+ 1 ) *DCOS ( -X ) / ( NHARD+ 1 ) 
Sc  -DSIN(-X)*(SIG1(1)*(U1(2)+DU1(1)  ) 

Sc  +SIG1(3)*(U1(1)+DU1(2))) 

Sc  -DCOS  ( -X)  *  (NHARD*EIGEN+1)  *  (SIGl  ( 1 )  *U1  ( 1 )  +SIG1  ( 3 )  *U1  (2 )  ) 


10  F0RMAT(1X,5(1X,E14.7) ) 

RETURN 

END 

************************************************************** 

SUBROUTINE  SAVE(X,yiNC) 


IMPLICIT 
INTEGER 
REAL*8 
REAL* 8 
REAL* 8 


NONE 

I 

X 

Y1NC(5) ,y(5) ,DyDX(5) ,MIS11(5) 
SIG1(4) ,DSIG1(4) ,DDSIG1(4) 


INTEGER  NHARD, COUNT 

REAL*8  RLAMBDA, SI, ALPHA. BETA, MAXMIS 

REAL*8  POTENTIALl (7, 2000) .CONSTANT! (3,4000) 

COMMON/MAINC/Sl , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/ POTENTIAL! . CONSTANT! 

COUNT=COUNT+! 

DO  I=!,5,l 

Y(I)=YINC(I) 

ENDDO 

POTENTIAL! ( ! , COUNT ) =X 
POTENTIAL! (7 , COUNT) =YINC ( 5 ) 

CALL  FCN  (5,X,Y,DYDX) 

POTENTIAL! ( 6 , COUNT ) =DYDX ( 4 ) 

DO  I=!,4,! 

POTENTIAL! ( I+! , COUNT) =Y ( I ) 

ENDDO 

Y ( 5 ) =DYDX ( 4 ) 

CALL  STRESS  (YINC (5 ) , Y. SIG!) 

CALL  DSTRESS  (YINC { 5 ) , Y, DSIG! ) 

CALL  DDSTRESS  (YINC (5) , Y, DDSIG!) 

CALL  SIGMIS (SIG! , SIG! , DSIG! , DSIG! , DDSIG! , DDSIG! , MIS!! ) 
MIS!!(!)=MIS!!(!)**0.5 
IF(MAXMIS.LT,MIS!!(!) )THEN 
MAXMIS=MIS!!(!) 

ENDIF 

FORMAT  (!X,6E!2.4) 

RETURN 

END 

SUBROUTINE  OUTPUT 

************************************************************** 


IMPLICIT 

INTEGER 

REAL*8 

REAL*8 

REAL*8 

REAL*8 


NONE 

I, J.KOUNT 
X.MEAN! 

YINC(5) ,Y(5) ,MIS!!(5) 

SIG!(4) ,DSIG!(4),DDSIG!(4) ,SS!(4) ,DSS!(4) ,DDSS!(4) 
EPS!(4) .DEPS!(4) ,DDEPS!(4) ,U!(2) ,DU!(2) 


INTEGER  NHARD , COUNT , lUl , IU2 , IU3 , IU4 , IU5 , IU6 , IU7 

REAL* 8  RLAMBDA, SI, ALPHA, BETA, MAXMIS 

REAL*8  POTENTIAL!(7, 2000) ,CONSTANT!(3, 4000) 

COMMON/MAINC/S! , RLAMBDA, ALPHA, BETA, NHARD, COUNT, MAXMIS 
COMMON/MAINA/ POTENTIAL! , CONSTANT! 

COMMON/ OUT/ lU!, IU2, IU3, IU4, IU5, IU6, IU7 


KOUNT=! 

DO  1=!, COUNT,! 

X=POTENTIAL! ( ! ,  I ) 
DO  J=!,5,! 
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c 


c 

c 


c 


c 


c 

c 

c 

c 


c 

c 

c 


y (J) =P0TENTIAL1 ( J+1, 1) 

ENDDO 

CALL  STRESS  {POTENTIALl (7, 1) , Y,  SIGl ) 

CALL  DSTRESS  (POTENTIALl (7 , 1) ,  Y,  DSIGl ) 

CALL  DDSTRESS  (POTENTIALl (7 , I) , Y, DDSIGl ) 

CALL  SIGMIS (SIGl , SIGl . DSIGl , DSIGl , DDSIGl , DDSIGl , MISll ) 

MIS11(1)=:MIS11(1)**0.5 

CALL  DEV (SIGl, SSI, MEANl) 

CALL  DDEV( DSIGl, DSSl) 

CALL  DDDEV( DDSIGl, DDSSl) 

CALL  STRAIN (MISll, SSI, EPSl) 

CALL ■ DSTRAIN (MISll , SSI , DSSl ,  DEPSl ) 

CALL  DDSTRAIN (MISll , SSI , DSSl , DDSSl , DDEPSl ) 

CALL  DISPLACE(POTENTIALl(7, I) , MISll, EPSl, DEPSl, Ul) 

CALL  DDISPLACE( POTENTIALl (7, 1)  ,MIS11 ,  EPSl , DEPSl , DDEPSl , DUl ) 
CALL  INTEGRAND ( X , POTENTIALl ( 7 , 1 ) ,  MIS 1 1 ,  S IGl , Ul , DUl , KOUNT) 
KOUNT=KOUNT+2 


WRITE (IU2, 15) 
WRITE (IU3, 10) 
WRITE (IU4, 10) 
WRITE (IU5, 10) 
WRITE (IU6, 10) 


X, (y(J),J=l, 5,1), POTENTIALl (7, I) 
X, (SIGl (J),J=1, 4,1) ,MIS11(1) 

X, (SSI (J),J=1, 4,1), MEANl 
X, (EPS1(J),J=1,4,1) 

X,U1(1) ,U1(2) ,DU1(1) ,DU1(2) 


ENDDO 


WRITEdl,  *) 
WRITE ( 11, *) 
WRITEdl,  *) 
WRITEdl,  *) 
WRITEdl,  *) 
WRITEdl,*) 


‘For  faster  convergence  use  the  following  guesses;' 


•  GUESS (1) 
■GUESS (2) 
•GUESS (3) 
•GUESS (4) 
■GUESS(5) 


• , POTENTIALl (4,1)/ POTENTIALl (2,1) 

', POTENTIALl (7,1) 

• , POTENTIALl ( 4 , COUNT) /POTENTIALl (2,1) 
• , POTENTIALl ( 5 , COUNT) /POTENTIALl (2,1) 
• , POTENTIALl (7 , COUNT) 


10  FORMAT  (lX,6E18dO) 
15  FORMAT  (lX,7E18dO) 


RETURN 

END 


SUBROUTINE  SORT ( NREC , ARRAY , NDIMl , NDIM2 ) 

************************************************************** 


IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  TEMP (6) 

DIMENSION  ARRAY (NDIMl, NDIM2) 

Sort  data 

DO  1=2, NREC, 1 
DO  J=2,NREC,1 

IF(ARRAY(1,  J)  .LT.ARRAYd,  J-1)  )THEN 
DO  K=l, NDIMl,  1 

TEMP ( K )= ARRAY ( K , J - 1 ) 

ENDDO 


oooo  on  n  non  o  noon 
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C 


DO  K=1,NDIM1,1 

ARRAY (K, J-1 ) =ARRAY (K, J ) 

ENDDO 

DO  K=1,NDIM1.1 

ARRAY (K,J)=TEMP(K) 

ENDDO 

.ENDIF 

ENDDO 

ENDDO 

RETURN 

END 

■k-k************************************************************ 

SUBROUTINE  INTEGCONST { COUNT , CONSTANTl ) 

*********•******•*************«*•«****«•****'***********«****** 


IMPLICIT  NONE 

INTEGER  I , COUNT 

REAL*8  CONSTANTl (3. 4000) 


INTEGER  lUl , IU2 , IU3 , IU4 , IU5 , IU6 , IU7 
COMMON/OUT/ lUl , IU2 , IU3 , IU4 . IU5 , IU6 , IU7 


Integrate  data 

DO  1=1, COUNT, 1 

CONSTANTl ( 3, 1)=0 
ENDDO 

DO  1=2, COUNT, 1 

CONSTANTl (3,1)= (CONSTANTl (2,1-1) +CONSTANT1 (2 , I ) ) /2 
Sc  *  (CONSTANTl  (1,1) -CONSTANTl  ( 1, 1-l )  ) 

Sc  +CONSTANTl(3, 1-l) 

ENDDO 

WRITE(IU7,10)  CONSTANTl ( 1, 1) , CONSTANTl (2,1) , 

Sc  CONSTANTl  (3, 1) ,  CONSTANTl  (3,  COUNT) 

DO  1=2, COUNT, 1 

WRITE (IU7, 10)  CONSTANTl (1,1) .CONSTANTl (2,1) , CONSTANTl ( 3 , 1 ) 
ENDDO 

10  F0RMAT(1X,4E15.7) 

RETURN 

END 

SUBROUTINE  OPENUNIT ( ITERATION, lUl , IU2 , IU3 , IU4 , IU5 , IU6 , IU7 ) 

*1i*1t*****ic1r************************ii***********l,ir***1i******i,ir1c 


IMPLICIT  NONE 

INTEGER  ITERATION, lUl , IU2 , IU3 , IU4 ,  IU5 , IU6 , IU7 
C 

lUl  =  11+ (ITERATION-1) *7  !  Program  output  file 
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c 


c 


c 


c 


c 

c 

c 


c 


IU2  =  12+ (ITERATION-1) *7 
IU3  =  13+ (ITERATION-1) *7 
IU4  =  14+ ( ITERATION- 1) *7 
IU5  =  15+ (ITERATION-1) *7 
IU6  =  16+ (ITERATION-1) *7 
IU7  =  17+ ( ITERATION- 1) *7 


Potential  and  derivatives 
Stresses 

Deviatoric  and  mean  stresses 
Strains 

displacements  and  derivatives 
integrand  and  constant 


OPEN  (UNIT=1U1 , STATUS= • NEW ‘ ) 
OPEN  (UNIT=IU2,STATUS=‘NEW ) 
OPEN  (UNIT=IU3,STATUS=‘NEW ) 
OPEN  (UN1T=1U4,STATUS=‘NEW ) 
OPEN  (UNIT=IU5,STATUS=*NEW ) 
OPEN  {UNIT=1U6, STATUS= 'NEW ) 
OPEN  {UNIT=IU7,STATUS=‘NEW ) 


WRITE (IU1,1) 

WRITE ( IU2 , 2 ) 

WRITE ( IU3 , 3 ) 

WRITE(IU4, 4) 

WRITE(IU5,5) 

WRITE(IU6, 6) 

WRITE ( IU7 , 7 ) 

1  FORMAT ( IX ' Program  output  * ) 

2  FORMAT(5X, 'Theta', 7X, 'fO',9X,  'fl'.9X,  '  f 2 ' , 

&  9X, ' f 3 ' , 9X, ' f 4 ' , 3X, ' eigenvalue ' ) 

3  FORMAT (5X, 'Theta' , 8X, 'sigrr ' , 7X,  ' sigtt ' , 7X, ' sigrt ' , 

St  7X, 'sigzz' ,  6X, 'sigmis' ) 

4  FORMAT (5X, 'Theta' , 6X, 'devsigrr ' , 4X, 'devsigtt ' , 

Sc  4X,  '  devsigrt ' ,  4X,  'devsigzz' ,  5X,  'meansig' ) 

5  FORMAT ( 5X, ' Theta ' , 7X, ' epsrr ' , 7X, ' epstt ' , 

&  7X, 'epsrt ' , 7X, 'epszz ' ) 

6  FORMAT(5X, 'Theta' ,8X, 'Ur'.lOX, 'Ut' ,10X, 'dUr' ,9X, 'dUt' ) 

7  FORMAT(7X, 'Theta', 8X, ' integrand' , 6X, 'integral', 

St  7X,  '  constant ' ) 


RETURN 

END 

*********ir*1,'tr******1r***********************1r*****1r**mr**ir***** 

SUBROUTINE  NORMALIZE (COUNT, MAXMIS, POTENTIALl ) 

************************************************************** 


IMPLICIT  NONE 
INTEGER  I, J, COUNT 
REAL*8  MAXMIS 
REAL*8  POTENTIALl (7, 2000) 

DO  1=1, COUNT, 1 
DO  J=2,6,l 

POTENTIALl (J, I) =POTENTIALl (J, I) /MAXMIS 
ENDDO 
ENDDO 


C 


RETURN 

END 


o  o  o  o 
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************************************************************** 

SUBROUTINE  CLOSEUNIT ( lUl , IU2 , IU3 , IU4 , IU5 , IU6 , IU7 ) 

************************************************************** 


IMPLICIT  NONE 

INTEGER  lUl, IU2, IU3 , IU4, IU5, IU6, IU7 
C 

CLOSE (UN1T=IU1) 

CLOSE (UNIT=IU2) 

CLOSE (UNIT=IU3) 

CLOSE (UNIT=IU4) 

CLOSE (UNIT=IU5) 

CLOSE (UNIT=IU6) 

CLOSE (UNIT=IU7) 

C 

RETURN 

END 

C 
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APPENDIX  D 

SHAPE  FUNCTIONS  USED  FOR  FINITE  ELEMENT  SOLUTIONS 
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Shape  functions  for  eight  node  quadratic  element; 

“  4  ^2) 

V2fe.«2)  =  :j(l  +  «j)(l-«2)(-l  +  «;-«2) 

Vj(li>?2)  “  4!*  ■'■  ■*■  ^2)(“'  ^2  ■'■  ^2) 

❖.(f;.«2)  =  :j(l-«;)(>  +  «2)H-«2  +  «2) 

V5(«,.l2)  =  |(l-«f){>-|2) 

'i'6(«2.«2)  =  |(l-«i)(l-«2) 


Coordinates  of  3x3  Gauss  quadrature  integration  points  and  corresponding 
quadrature  weights: 


(l/.4|)=(o.o).  wj=|j 


APPENDIX  E 


AB  AQUS  FINITE  ELEMENT  INPUT  FILE 
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♦HEADING 

CCT_AW_01_3 00000 

**  NEUTRAL  FILE  GENERATED  ON:06-Jan-94  14:51:44  PATABA  VERSION:  3 . lA 

** 

**  NODE  DEFINITIONS 

** 

♦NODE 

1,  O.OOOOOOOOOE-00,  O.OOOOOOOOOE-00,  O.OOOOOOOOOE+00 

2,  -0.122028234E-11,  0 .591932803E-05,  0 . OOOOOOOOOE+00 

3,  O.OOOOOOOOOE-00,  O.OOOOOOOOOE-00,  0 . OOOOOOOOOE+00 


3688,  0.135783806E+00,  -0.135783881E+00,  0 . OOOOOOOOOE+00 

3689,  0.132622838E+00,  -0 . 161601558E+00 ,  0 . OOOOOOOOOE+00 

3690,  0.106684595E+00,  -0 . 159664989E+00 ,  0 . OOOOOOOOOE+00 

** 

♦♦  NODE  SETS  FROM  NAMED  COMPONENTS 
♦NSET,  NSET=CRACK 

1  3  5  6  8  9  11  12  14  15  17  18  20  21  23  24 

26  27  29  30  32  33  35  36  38  39  41  42  44  45  48  49 

50 


♦♦  ELEMENT  DEFINITIONS 

** 

♦ELEMENT,  TyPE=CPE8  ,  ELSET=PID1 

1,  2565,  2497,  2500,  2568,  2496,  2499,  2498,  2566 

2,  2622,  2565,  2568,  2625,  2567,  2566,  2569,  2623 

3,  2667,  2622,  2625,  2670,  2624,  2623,  2626,  2668 


1132, 

3638, 

3644, 

3605, 

3599, 

3640, 

3607,  3601,  3600 

1133, 

3645, 

3641, 

3602, 

3606, 

3642, 

3604,  3603,  3609 

1134, 

3644, 

3645, 

3606, 

3605, 

3647, 

3609,  3608,  3607 

♦ELEMENT,  TyPE=CPE8 

ELSET: 

=PID2 

865, 

2772, 

2770, 

2736, 

2738, 

2769, 

2735,  2734,  2737 

866, 

2775, 

2772, 

2738, 

2741, 

2773, 

2737,  2739,  2740 

867, 

2778, 

2775, 

2741, 

2744, 

2776, 

2740,  2742,  2743 

868, 

2781, 

2778, 

2744, 

2747, 

2779, 

2743,  2745,  2746 

869, 

2784, 

2781, 

2747, 

2750, 

2782, 

2746,  2748,  2749 

870, 

2787, 

2784, 

2750, 

2753, 

2785, 

2749,  2751,  2752 

871, 

2790, 

2787, 

2753, 

2756, 

2788, 

2752,  2754,  2755 

872, 

2796, 

2790, 

2756, 

2762, 

2791, 

2755,  2757,  2761 

873, 

2794, 

2799, 

2766, 

2760, 

2793, 

2764,  2759,  2758 

874, 

2800, 

2796, 

2762, 

2767, 

2797, 

2761,  2763,  2765 

875, 

2799, 

2800, 

2767, 

2766, 

2802, 

2765,  2768,  2764 

** 


*  * 


*  * 


ELEMENT  PROPERTIES 
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♦SOLID  SECTION, ELSET=PID1 , MATERIAL=PLASTIC,  ORIENTATION=CyLINDER 

1.0 

♦SOLID  SECTION, ELSET=PID2 , MATERIAL=ELASTIC, ORIENTATION=CYLINDER 

1.0 

** 

♦♦  DEFINE  A  CYLINDRICAL  COORDINATE  SYSTEM  TO  GET  OUTPUT 

** 

♦ORIENTATION, NAME=CYLINDER, DEFINITION=COORDINATES , SYSTEM=CYLINDRICAL 
0.0, 0.0, 0.0, 0.0, 0.0, 1.0 
** 

♦♦  MATERIAL  DEFINITIONS 
** 

♦MATERIAL, NAME=ELASTIC 

♦ELASTIC 

30000000,0.3 

** 

♦MATERIAL , NAME=PLASTIC 
♦USER  MATERIAL, CONSTANTS=4 
60000,0.002,10,0.3 


★  * 

♦♦  LOAD  CASE  1 

*  * 

♦STEP,AMPL1TUDE=RAMP, INC=500 

LOAD  CASE  1 

♦STATIC 

0.1,  1.0,  0.0000001,  1 

** 

♦DLOAD,  OP=NEW 

865,  P3,  -3750.0 

866,  P3,  -3750.0 

867,  P3,  -3750.0 

868,  P3,  -3750.0 

869,  P3,  -3750.0 

870,  P3,  -3750.0 

871,  P3,  -3750.0 

872,  P3,  -3750.0 

873,  P3,  -3750.0 

874,  P3,  -3750.0 

875,  P3,  -3750.0 

♦♦  Zero  displacements  on  planes  of  symmetry 

it  * 

♦BOUNDARY,  OP=NEW 

10,  1,,  0.0 

12,  1,,  0.0 

61,  1,,  0.0 
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3615,  1,,  0.0 

3617,  1,,  0.0 

3618,  1,,  0.0 

**  Fracture  Mechanics 

★  ★ 

* J-INTEGRAL, CONTOURS=50 , FREQUENCY=1 , OUTPUT=BOTH, SYMM 

1.0, 0.0 

CRACK 

** 

**  Output 

•PRINT, CONTACT=NO, FREQUENCY=1 ,  RESIDUAL=YES , SOLVE=YES 

** 

♦NODE  FILE,FREQUENCY=500,GLOBAL=YES 
U 

•NODE  PRINT, FREQUENCY=500,GLOBAL=YES 

U 

RF 

•EL  FILE,  FREQUENCY=500,POSITION=INTEGRATION  POINTS 
S 

SINV 

E 

•EL  PRINT, FREQ= 5 00,POS= INTEGRATION  POINTS,  SUMMARY=NO, TOTALS=NO 
S 

SINV 

E 

•FILE  FORMAT,  ASCII 
•END  STEP 

** 


•STEP, AMPLITUDE=RAMP, INC=500 

LOAD  CASE  2 

♦STATIC 

0.1,  1.0,  0.0000001,  1 


♦DLOAD,  OP=NEW 

865,  P3,  -7500.0 

866,  P3,  -7500.0 

867,  P3,  -7500.0 

868,  P3,  -7500.0 

869,  P3,  -7500.0 

870,  P3,  -7500.0 

871,  P3,  -7500.0 

872,  P3,  -7500.0 

873,  P3,  -7500.0 

874,  P3,  -7500,0 

875,  P3,  -7500.0 

** 

••  Fracture  Mechanics 


•J- INTEGRAL, CONTOURS=50 , FREQUENCY=1 , OUTPUT=BOTH , SYMM 
1.0, 0.0 
CRACK 


*  * 


*  *  Output 
** 

•PRINT, CONTACT=NO. FREQUENCY =1 , RESIDUAL=yES, SOLVE=YES 

** 

•NODE  FILE,FREQUENCY=500,GLOBAL=YES 
U 

•NODE  PRINT,FREQUENCY=500,GLOBAL=YES 

U 

RF 

•EL  FILE,  FREQUENCY=500,POSITION=INTEGRATION  POINTS 
S 

SINV 

E 

•EL  PRINT, FREQ=500, POS=INTEGRATION  POINTS, SUMMARY =NO, TOTALSaNO 
S 

SINV 

E 

•FILE  FORMAT,  ASCII 
•END  STEP 


•STEP,AMPLITUDE=RAMP, INC=500 

LOAD  CASE  20 

•STATIC 

0.1,  1.0,  0.0000001,  1 

** 


UO,  OP=NEW 

865, 

P3, 

-75000.0 

866, 

P3, 

-75000.0 

867, 

P3, 

-75000.0 

868, 

P3, 

-75000.0 

869, 

P3, 

-75000.0 

870, 

P3, 

-75000.0 

871, 

P3, 

-75000.0 

872, 

P3, 

-75000.0 

873, 

P3, 

-75000.0 

874, 

P3, 

-75000.0 

875, 

P3, 

-75000.0 

••  Fracture  Mechanics 


•J- INTEGRAL, CONTOURS=50 , FREQUENCY=1,  OUTPUT=BOTH, SYMM 

1.0, 0.0 

CRACK 


★  ★ 


Output 


‘PRINT, CONTACT=NO, FREQUENCY=1 , RESIDUAL=YES , SOLVE=YES 

** 

♦NODE  FILE, FREQUENCY=500 , GLOBAL=YES 
U 

‘NODE  PRINT, FREQUENCY=500,GLOBAL=YES 

U 

RF 

*EL  FILE,  FREQUENCY=500,POSITION=INTEGRATION  POINTS 
S 

SINV 

E 

*EL  PRINT, FREQ=500,POS=INTEGRATION  POINTS,  SUMMARY=NO, TOTALS s:NO 
S 

SINV 

E 

‘FILE  FORMAT,  ASCII 
‘END  STEP 


APPENDIX  F 


INTEGRATION  CONSTANTS  FOR  ANALYTICAL  /-INTERGAL 
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APPENDIX  G 


FORTRAN  PROGRAM  TO  EXTRACT  AND  NORMALIZE  THE 
FINITE  ELEMENT  STRESSES  AND  DETERMINE 


THE  ANALYTICAL  AMPLITUDES 


onn  non  noon  nnon  non  ooo  ooo  non  n  n  n  on 
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PROGRAM  EXTRACTSTRESS 

*************************************************************** 


IMPLICIT  REAL*8  (A-H,0-Z) 

CHARACTER*80  FILE ( 50 ) , FIELDFILE ( 10 ) 

DIMENS ION  OUTPUTl { 4 ) , OUTPUT2 ( 4 ) , RJINT (100), NUMANGLES (10) 
DIMENSION  EIGEN(IO) . ATWO (3 ) . AFOUR (3 ) , SDIFFl (3 ) 

DIMENSION  NRAYS (17,300), COORDINATES (5000,4), VARIABLES (5000,100) 
DIMENSION  NRINGS(55,300) ,ASYMSTRESS(10, 10) ,RL(10, 3) ,RG(6) 
DIMENSION  STRESSFIELD(10,2000,10) 

Enter  needed  constants 

SIGNOT=60000 
EM0D=3 0000000 
EPSNOT=0.002 

Get  the  names  of  the  PATRAN.NOD  files  with  stresses 
CALL  FILENAMES (NFILES , FILE) 

Get  the  groups  of  nodes  for  the  individual  rays 
CALL  RAYS (NSET1,NUMN0DES1, NRAYS) 

Get  the  groups  of  nodes  for  the  individual  rings 
CALL  RINGS (NSET2 , NUMNODES2 , NRINGS ) 

Get  the  coordinates  of  the  nodes 
CALL  COORDS (NUM, COORDINATES) 

Get  the  J-integral  values  for  the  load  steps  corresponding 
to  the  results  in  the  PATRAN.NOD  files 

CALL  JINTEGRALS(NUMJINT, RJINT) 

Get  the  filenames  which  the  analytic  angular  solutions 
are  stored  in 

CALL  FIELDFILENAMES (NFIELDFILES, FIELDFILE) 

Read  in  the  analytic  angular  solutions 

CALL  FIELDS (NFIELDFILES , FIELDFILE, NUMANGLES , STRESSFIELD) 

Get  the  eigenvalues  for  the  analytic  solution 
CALL  EIGENVALUES (NEIGEN, EIGEN,  FIRSTCONSTANT) 

DO  J=l, NFILES,! 

CALL  VARIABLE(FILE(J) , NUM, VARIABLES) 

DO  MMM=1,3,1 


C 


o  o  o 
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WRITE (20, *) 

WRITEOO,  *) 

ENDDO 

WRITE(20,2)FILE(J) 

WRITE(30,2)FILE(J) 

Find  variables  for  given  nodes  in  node  sets  in  rays 


& 


& 

Sc 


Sc 


Sc 

Sc 

Sc 

Sc 

Sc 

Sc 


DO  K=1,NSET1,1 
WRITE(20,  *) 

WRITE (20, *) 'Ray:  '.K 
WRITE(20, *) 

WRITE(20,*) 'Nonnalized_Quantities : ' 

WRITE(20, *) 

WRITE (20, 3) 

WRITE(20,  *) 

WRITE(30,  *) 

WRITE(30, *) 'Ray: '  ,K 
WRITE (30, *) 

WRITE(30,*) ‘Normali2ed_Quantities : ' 

WRITE(30, *) 

WRITE(30,4) 

WRITE(30, *) 

A2=0.0 

A4=0.0 

DO  I=l,NUMNODESl,l 

CALL  SEARCHRAYS ( I , NUM, K , COORDINATES , NRAYS , VARIABLES , 
RR, THETA, OUTPUTl ) 
RNORM=RR*SIGNOT/RJINT ( J ) 

IF( (RNORM.lt. 15.0) .AND. (RNORM.GT. 0 . 0 ) ) THEN 
IF(THETA.LT.3 .1415)THEN 

CALL  INTERPOLATESTRESS ( THETA , NFIELDFILES , 

NUMANGLES, STRESSFIELD, ASYMSTRESS) 
CALL  COEFFl (SIGNOT,  EPSNOT, FIRSTCONSTANT, 

ASYMSTRESS, RR,RJINT(J) , EIGEN, RL) 
CALL  COEFF2 (RL, RG, OUTPUTl ( 1 ) , OUTPUTl (2 ) ) 

CALL  NEWTON (RG, RL, OUTPUTl ( 1 ) , OUTPUTl ( 4 ) , 

A2 , A4 , ATWO , AFOUR , SDIFFl ) 

ENDIF 

RLAMBDAsOUTPUTl ( 3 ) / ( OUTPUTl ( 1 ) +OUTPUT1 ( 2 ) ) 
WRITE(20, 1)  RR*SIGNOT/RJINT(J) , THETA, 

OUTPUTl ( 1 ) /SIGNOT, 

OUTPUTl (2) /SIGNOT, 

OUTPUTl (3 ) /SIGNOT, 

OUTPUTl ( 4 ) /SIGNOT, RLAMBDA 
WRITE (30, 5)  RR*SIGNOT/RJINT(J) , THETA, RLAMBDA, 
(ATWO(LL) .LL=1,3, 1) , 

(AFOUR(LL) ,LL=1,3, 1) 

ENDIF 

ENDDO 

ENDDO 


Find  variables  for  given  nodes  in  node  sets  in  rings 


C 

C 

C 


DO  NNN=1,2, 1 
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Sc 


Sc 

Sc 


Sc 


Sc 

Sc 

Sc 

Sc 


WRITE (20, *) 

WRITEOO,  *) 

ENDDO 

WRITE{20, *) 

WRITE(20,*) 'Ring;  First_with_nonnalized_radius_less_than_15 ' 
WRITE(20,*) 

WRITE{20, *) 'Nonnali2ed_Quantities : ' 

WRITE(20, *) 

WRITE(20,3) 

WRITE(20,*) 

WRITEOO,*) 

WRITEOO,*)  'Ring:  First_with_nonnali2ed_radius_less_than_15 ' 
WRITEOO,  *) 

WRITEOO,*)  'Nonnali2ed_Quantities :  ' 

WRITEOO,  *) 

WRITEOO, 4) 

WRITEOO,*) 

DO  K=1,NSET2,1 
A2=0.0 
A4=0.0 

DO  I=l,NUMNODES2,l 

CALL  SEARCHRINGS ( I ,  NUM,  K, COORDINATES , NRINGS , 

VARIABLES , RR , THETA , OUTPUT2 ) 
RNORM=RR*SIGNOT/RJINT { J ) 

IF{ {RNORM.LT.15.0) .AND. (RNORM.GT. 0 . 0 ) ) THEN 
IF(THETA.LT.3.1415)THEN 

CALL  INTERPOLATESTRESS ( THETA , NFIELDFILES , 

NUMANGLES , STRESSFIELD, ASYMSTRESS ) 
CALL  COEFFl (SIGNOT,  EPSNOT, FIRSTCONSTANT, 

ASYMSTRESS , RR , RJ INT ( J ), EIGEN , RL ) 
CALL  COEFF2 ( RL , RG , OUTPUT2 ( 1 ) , OUTPUT2 ( 2 ) ) 

CALL  NEWTON ( RG , RL , OUTPUT2 ( 1 ) , OUTPUT2 ( 4 ) , 

A2 , A4 , ATWO , AFOUR , SDIFFl ) 

ENDIF 

RLAMBDA=OUTPUT2 ( 3 ) / ( OUTPUT2 ( 1 ) +OUTPUT2 ( 2 ) ) 

WRITE (2 0,1)  RR*SIGNOT/RJINT(J) , THETA, 

OUTPUT2 ( 1 ) / S IGNOT , 

OUTPUT2 ( 2 ) /SIGNOT, 

OUTPUT2 (3 ) /SIGNOT, 

OUTPUT2 ( 4 ) /SIGNOT, RLAMBDA 
WRITE (3 0,5)  RR*SIGNOT/RJINT(J) , THETA, RLAMBDA, 

( ATWO (LL),LL=1, 3,1)  , 

(AFOUR {LL),LL=1, 3,1) 

ENDIF 

ENDDO 

IF ( (RNORM.lt. 15.0) .AND. (RNORM.GT. 0 . 0 ) )THEN 
IF(K.LE. (NSET2-2) )THEN 
WRITE(20,*) 

WRITE(20, *) 'Ring: ' ,K+1 
WRITE(20,*) 

WRITE(20,*) 'Nonnali2ed_Quantities ; ' 

WRITE(20,*) 

WRITE (20, 3) 

WRITE (20,*) 

WRITE(30, *) 


onoo  oof)  o  no  oono 
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WRITEOO,  *)  'Ring;  ‘  ,K+1 
WRITEOO,  *) 

WRITE{30, *) •Nonnalized_Quantities: ' 

WRITEOO,  *) 

WRITEOO, 4) 

WRITEOO,  *) 

ENDIF 
ENDIF 
ENDDO 

ENDDO 

1  FORMAT(10(TR1,E12.4) ) 

2  FORMAT (A80) 

3  F0RMAT(TR6, 'rbar' ,TR8, 'theta*  ,TR8,  'sigrr' ,TR8, 'sigtf ,TR8, 

&  'sigzz' ,TR8, 'sigrt' ,TR7, 'lambda' ) 

4  FORMAT {TR6 , ' rbar ' , TR8 , ' theta ' ,  TR7 , ' lambda ' , TR9 , ' A2_l ' , TR9 , 

&  ‘A2_2',TR9, 'A2_3',TR9, 'A4_1',TR9, 'A4_2',TR9, ‘A4_3') 

5  F0RMAT(9(TR1,E12.4) ) 

STOP 
END 

★  ♦★★**★*★♦★*★♦*♦★★★**•★★★♦★*♦♦♦*★★*★★★***★**♦♦*★★★*★*♦♦**★★*★★★* 

SUBROUTINE  FILENAMES (NFILES, FILE) 

★  ★♦★*★♦★***♦'★★★★★**★*★**♦★*★***★*★♦★*'>►•*★★*******★*★★★★'*♦★★*★★★ 

IMPLICIT  REAL*8  (A-H,0-Z) 

CHARACTER* 80  FILE (50) 

OPEN (UNIT=14 , FILE= ' filenames . dat ' ,  STATUS= ' old ' ) 

KKK=1 

100  READ(14,1,END=200)FILE(KKK) 

KKK=KKK+1 
GOTO  100 

200  KKK=KKK-1 
NFILES=KKK 

CLOSE (UNIT=14) 

1  FORMAT (A80) 

RETURN 
END 

SUBROUTINE  FIELDFILENAMES (NFIELDFILES, FIELDFILE) 

IMPLICIT  REAL*8  {A-H,0-Z) 

CHARACTER*80  FIELDFILE ( 10 ) 

C 

OPEN (UNIT=16 , FILE= ' f ieldf ilenames . dat  * ,  STATUS= ' old ' ) 
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c 

KKK=1 

100  READ(16,1,END=200)FIELDFILE(KKK) 

KKK=KKK+1 
GOTO  100 
C 

200  KKK=KEK-1 

NFIELDFILES=KKK 

C 

CLOSE (UNIT=16) 

C 

1  FORMAT (A80) 

C 

RETURN 

END 

C 

^  ********•****'****************************•*****************«** 

SUBROUTINE  RAYS  (NSETl, NUl'iNODESl,  NRAYS) 

^  ************•**•********•**»****««********«'«*****•**•***«***** 

C 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  NRAYS(17,300) 

C 

OPEN {UNIT=13 , FILE= ' rays . dat • , STATUS= ' old ' ) 

C 

1  FORMAT(16(I5) ) 

C 

READ ( 13 , * ) NUMNODESl 
III»1 

100  READ(13,1,ERR=200,END=300) (NRAYS ( III, J) , J=1 , NUMNODESl ) 
III=III+1 
200  GOTO  100 
300  III=III-1 
NSET1=III 
C 

CLOSE (UNIT=13) 

C 

RETURN 

END 

C 

Q  ***************************************************************** 

SUBROUTINE  RINGS ( NSET2 , NUMNODES2 ,  NRINGS ) 

C  ************************************************************** 

c 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  NRINGS(55,300) 

C 

OPEN {UNIT=13 , FILE= ’ rings . dat ' , STATUS= ' old ’ ) 

C 

1  FORMAT (16 (15) ) 

C 

READ { 13 , * ) NUMNODES2 
III=1 

100  READ(13, 1,ERR=200,END=300) (NRINGS ( III , J) , J=1 , NUMNODES2 ) 
III=III+1 


200  GOTO  100 
300 

NSET2=III 
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c 

CLOSE {UNIT=13) 

C 

RETURN 

END 

C 

Q  A***'**************************************************'******** 

SUBROUTINE  COORDS ( NUM , COORDINATES ) 

0  ************************************************************** 

C 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  COORDINATES (5000, 4) 

C 

OPEN (UNIT=10 , FILE= ■ coordinates . daf  ,  STATUS= ■ old ' ) 

C 

1  FORMAT ( 18 , IX, E17 . 9 , IX, E17 . 9 , IX, E17 . 9 ) 

C 

NUM=1 

100  READdO,  1,ERR=200)  NODENUMBER, COORDINATES  (NUM,  2 )  , 

&  COORDINATES (NUM, 3 ) , COORDINATES (NUM, 4 ) 

COORDINATES (NUM, 1 ) =NODENUMBER 
NUM=NUM+1 
GOTO  100 
200  NUM=NUM-1 
C 

CLOSE (UNIT=10) 

C 

RETURN 

END 

C 

c  *************************************************************** 

SUBROUTINE  VARIABLE (NAME, NUM, VARIABLES) 

Q  *********1t*****************1t**1r*****1r************1r1,1r*****tiii*** 

C 

IMPLICIT  REAL*8  (A-H,0-Z) 

CHARACTER* 80  NAME 
DIMENSION  VARIABLES (5000, 100) 

C 

1  FORMAT(I8, 5(E13.7) ) 

2  FORMAT(5(E13.7) ) 

3  FORMAT (E13. 7) 

C 

OPEN (UNIT=11 , FILE=NAME, STATUS= • old • ) 

C 

READdl,  *) 

READdl,  *) 

READdl,  *) 

READdl,*) 

DO  1=1, NUM, 1 

READ (11,1)  NODENUMBER, (VARIABLES ( I , J) , J=2 ,6,1) 

VARIABLES (1,1) =NODENUMBER 
DO  K=l,17,l 
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NSTART=7+ (K-l) *5 
NEND=11+(K-1)*5 

READ (11,2)  (VARIABLES ( I , J ) , J=NSTART, NEND, 1 ) 

ENDDO 

READ (11, 3)  VARIABLES(I,NEND+1) 

ENDDO 

C 

CLOSE  ■(UNIT=11) 

C 

RETURN 

END 

C 

^  **********«********************•***************'*************** 

SUBROUTINE  J INTEGRALS (NUMJINT, RJINT) 

Q  ************************************************************** 

c 

c 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  RJINT (100) 

C 

OPEN (UNIT=15,FILE=‘j integral. daf ,STATUS=‘ old* ) 

C 

III=1 

100  READ(15, 1,END=200)  RJINT(III) 

III=III+1 
GOTO  100 
C 

200  III=III-1 
NUMJINT=III 
C 

CLOSE (UNIT=15) 

C 

1  FORMAT (FIO. 2) 

C 

RETURN 

END 

C 

C  ************************************************************** 

SUBROUTINE  SEARCHRAYS ( I , NUM, NSETl , COORDINATES , NRAYS , 

&  VARIABLES , RR , THETA, OUTPUTl ) 

C  ************************************************************** 

c 

c 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  OUTPUTl (4) 

DIMENSION  COORDINATES (5000, 4) ,NRAYS(17, 300) , VARIABLES (5000, 100) 
C 

C  Find  node  coordinates 

C 

ZERO=0. 00000000000 
ONE=l. 000000000000 
PIE=4.0*ATAN(ONE) 

DO  J=1,NUM, 1 

IF (NRAYS (NSETl, I) .EQ. INT( COORDINATES (J,  1 ) ) )THEN 
XX=COORDINATES ( J, 2 ) 


oon  no  o  ooo 
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YY=C00RDINATES ( J, 3 ) 
RR={XX**2.0+YY**2.0)**0.5 
IF (YY.EQ. ZERO) THEN 
THETA=PIE 
GOTO  100 
ENDIF 

THETA=ATAN(ABS(XX)/ABS(YY) ) 

IF (YY .GT. ZERO) THETA= PIE-THETA 
GOTO  100 
ENDIF 
ENDDO 


Find  variables 
100  DO  J=1,NUM, 1 

IF(NRAYS(NSET1, I) .EQ. INT( VARIABLES (J, 1 ) ) )THEN 
OUTPUTl ( 1 ) =VARIABLES ( J,  2 ) 

OUTPUTl { 2 ) =VARIABLES ( J , 3 ) 

OUTPUTl (3 ) =VARIABLES ( J, 4) 

OUTPUTl ( 4 ) =VARIABLES { J,  5 ) 

GOTO  200 
ENDIF 
ENDDO 

200  RETURN 
END 


•****•*****•«***•*«***•*******«•*•********••«•**********«****• 

SUBROUTINE  SEARCHRINGS { I , NUM,  NSET2 , COORDINATES , NRINGS , 

&  VARIABLES , RR , THETA , OUTPUT2 ) 

*************************************************************** 


IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  OUTPUT2(4) 

DIMENSION  COORDINATES (5000,4), NRINGS (55,300) 
DIMENSION  VARIABLES (5000, 100) 

C 

C  Find  node  coordinates 
C 

ZERO=0. 00000000000 
ONE=l. 000000000000 
PIE=4.0*ATAN(ONE) 

DO  J=1,NUM,1 

IF ( NRINGS ( NSET2 , I ) . EQ . INT ( COORDINATES ( J , 1 ) ) ) THEN 
XX=COORDINATES ( J, 2 ) 

YY=COORDINATES ( J,  3 ) 

RR=(XX**2.0+YY**2.0)**0.5 
IF  ( YY .  EQ .  ZERO )  THEN 
THETA=PIE 
GOTO  100 
ENDIF 

THETA=ATAN(ABS(XX) /ABS(YY) ) 

IF ( YY .GT. ZERO) THETA=PIE-THETA 
GOTO  100 


ooo  oo  on  nonoo  n  ooooo  n  n  n  n 
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ENDIF 

ENDDO 


Find  variables 
100  DO  J=1,NUM,1 

IF(.NRINGS(NSET2,I)  .EQ. INT( VARIABLES  (J,  1 )  )  )THEN 
OUTPUT2 ( 1 ) =VARIABLES { J, 2 ) 

OUTPUT2 ( 2 ) =VARI ABLES { J , 3 ) 

OUTPUT2 ( 3 ) =VARI ABLES ( J , 4 ) 

OUTPUT2 ( 4 ) =VAR TABLES { J , 5 ) 

GOTO  200 
ENDIF 
ENDDO 

200  RETURN 
END 

*************************************************************** 

SUBROUTINE  FIELDS (NFIELDFILES, FIELDFILE, NUMANGLES, STRESSFIELD) 


IMPLICIT  REAL*8  (A-H,0-Z) 

CHARACTER* 80  FIELDFILE(IO) 

DIMENSION  NUMANGLES(IO) 

DIMENSION  STRESSFIELD{10,2000,10) 

DO  Irl,NFIELDFILES,l 

OPEN (UNIT=17 , FILE=FIELDFILE ( I ) ,  STATUSs ' old • ) 

KKK=1 

100  READ(17,1,ERR=200,END=300) (STRESSFIELD ( I , KKK, J) , J=1 , 6, 1 ) 

IF{I.EQ.6)THEN 
DO  LLL=2,6,1 

STRESSFIELD ( I , KKK, LLL) =0.0 
ENDDO 
ENDIF 
KKK=KKK+1 
200  GOTO  100 

300  KKK=KKK-1 

NUMANGLES(I)=KKK 
CLOSE (UNIT=17) 

ENDDO 

1  FORMAT{6(E18.10)  ) 

RETURN 

END 

************************************************************** 
SUBROUTINE  INTERPOLATESTRESS (THETA, NFIELDFILES , NUMANGLES , 

&  STRESSFIELD, ASYMSTRESS) 

♦************************************************************* 


OOCJCJ  UO  UCJO  u  u 
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IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  NUMANGLES(IO) 

DIMENSION  ASYMSTRESS(10,10) 

DIMENSION  STRESSFIELD(10,2000,10) 

C 

DO  I=1,NFIELDFILES,  1 

DO  J=2 , NUMANGLES ( I )  ,  1 

’lF( (STRESSFIELD(I,J-1,1) .LE. THETA) .AND. 

&  (STRESSFIELDd,  J.l)  .GE. THETA)  ) THEN 

DO  K=l,5,l 

ASYMSTRESS { I , K) =STRESSFIELD ( I , J-1 , K+ 1 ) + 

&  (THETA-STRESSFIELDd.  J-1,1)  )/ 

fic  ( STRESSFIELD (I , J, 1 ) -STRESSFIELD (I , J-1 , 1 ) ) * 

&  (STRESSFIELDd,  J,K+1 ) -STRESSFIELDd,  J-1, K+1)  ) 

ENDDO 
GOTO  100 
ENDIF 
ENDDO 
100  ENDDO 
C 

RETURN 
END 

SUBROUTINE  EIGENVALUES (NEIGEN, EIGEN, FIRSTCONSTANT) 

*•««**********«**•**•***'*«*******•*'*'****»*************•«*****• 

IMPLICIT  REAL*8  {A-H,0-Z) 

DIMENSION  EIGEN (10) 

OPEN (UNIT=18,FILE=' eigenvalues. daf  ,STATUS=‘ old' ) 

100  READ(18, 1,ERR=100,END=4 00) FIRSTCONSTANT 
KKK=1 

200  READ ( 18, 1,ERR=3 00, END=400) EIGEN (KKK) 

KKK=KKK+1 
300  GOTO  200 
400  KKK=KKK-1 
NEIGENxKKK 

CLOSE (UNIT=18) 

1  FORMAT (Fll. 7) 

RETURN 
END 

SUBROUTINE  COEFFl ( SIGNOT, EPSNOT, FIRSTCONSTANT, ASYMSTRESS ,  RR, 
t  RJINTEGRAL, EIGEN, RL) 

************************************************************** 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  EIGEN (10) 

DIMENSION  ASYMSTRESS ( 10, 10) ,RL( 10, 3) 
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c 

RBAR=SIGNOT*EPSNOT*RR/RJINTEGRAL 
BETA=RJINTEGRAL/ (SIGNOT*EPSNOT) 

C 

RL { 1 , 1 ) =SIGNOT*FIRSTCONSTANT**EIGEN ( 1 ) * ASYMSTRESS (1,1)* 

&  RBAR**EIGEN(1) 

RL ( 1 , ? ) =SIGNOT*FIRSTCONSTANT**EIGEN ( 1 ) *ASYMSTRESS (1,2)* 

&  RBAR**EIGEN(1) 

RL ( 1 , 3 ) =SIGNOT*FIRSTCONSTANT**EIGEN ( 1 ) *ASYMSTRESS (1,3)* 

Sc  RBAR**EIGEN(1) 

RL ( 2 , 1 ) =SIGNOT*BETA*  *EIGEN ( 2 ) *ASYMSTRESS (2,1) *RBAR*  *EIGEN ( 2 ) 

RL ( 2 , 2 ) =SIGNOT*BETA*  *EIGEN ( 2 ) * ASYMSTRESS (2,2) *RBAR*  *EIGEN ( 2 ) 

RL ( 2 , 3 ) =SIGNOT*BETA*  *EIGEN ( 2 ) * ASYMSTRESS (2,3) *RBAR*  *EIGEN ( 2 ) 
RL(3, 1) =SIGNOT*BETA**EIGEN(3) * (BETA/FIRSTCONSTANT) **EIGEN(1) * 

&  ASYMSTRESS (3, 1)*RBAR**EIGEN (3) 

RL(3 , 2) =SIGNOT*BETA**EIGEN(3) * (BETA/FIRSTCONSTANT) **EIGEN(1) * 

Sc  ASYMSTRESS  (3, 2)  *RBAR**EIGEN (3) 

RL ( 3 , 3 ) =SIGNOT*BETA*  *EIGEN ( 3 ) * (BETA/FIRSTCONSTANT) *  *EIGEN ( 1 ) * 

Sc  ASYMSTRESS  (3, 3  )*RBAR**EIGEN(3) 

RL  ( 4 , 1 ) =SIGNOT*BETA*  *EIGEN ( 4 ) * ASYMSTRESS (4,1) *RBAR*  *EIGEN ( 4 ) 

RL  ( 4 , 2 ) =S1GN0T*BETA**EIGEN ( 4 ) * ASYMSTRESS (4,2) *RBAR*  *EIGEN ( 4 ) 

RL ( 4 , 3 ) =SIGNOT*BETA**EIGEN ( 4 ) *ASYMSTRESS (4,3) *RBAR*  *EIGEN ( 4 ) 

RL { 5 , 1 ) =SIGNOT*BETA*  *EIGEN ( 5 ) * (BETA/FIRSTCONSTANT) 

Sc  *  *  ( 2  *EIGEN  ( 1 ) )  *  ASYMSTRESS  (5,1)  *RBAR*  *EIGEN  ( 5 ) 

RL ( 5 , 2 ) =SIGNOT*BETA*  *EIGEN ( 5 ) * ( BETA/FIRSTCONSTANT ) 

Sc  **(2*EIGEN(1)  )*ASYMSTRESS(5,2)*RBAR**EIGEN(5) 

RL(5, 3) =SIGNOT*BETA**EIGEN(5) * (BETA/FIRSTCONSTANT) 

Sc  **(2*EIGEN(1) )*ASYMSTRESS(5,3)*RBAR**EIGEN(5) 

RL(6, 1) =SIGNOT*BETA**EIGEN(6) * (BETA/FIRSTCONSTANT) **EIGEN(1) * 

Sc  ASYMSTRESS  (6, 1)*RBAR**EIGEN( 6) 

RL ( 6 , 2 ) =SIGNOT*BETA**EIGEN ( 6) * (BETA/FIRSTCONSTANT) *  *EIGEN ( 1 ) * 

&  ASYMSTRESS (6, 2 )*RBAR**EIGEN( 6) 

RL { 6 , 3 ) =SIGNOT*BETA*  *EIGEN ( 6 ) * ( BETA/FIRSTCONSTANT ) *  *EIGEN ( 1 ) * 

Sc  ASYMSTRESS  (6, 3  )*RBAR**EIGEN( 6) 

C 

RETURN 

END 

C 

Q  It************************************************************* 

SUBROUTINE  COEFF2 (RL, RG, SIGRR, SIGTT) 

^  ************************************************************** 

C 

IMPLICIT  REAL*8  {A-H,0-Z) 

DIMENSION  RG(6) 

DIMENSION  RL(10,3) 

C 

RG(1)=RL(6,2)*RL(5, 1)-RL(6,1)*RL(5,2) 
RG(2)=RL(4,2)*RL(5.1)+RL(6,2)*RL(3,1)-RL(4,1)*RL(5,2) 

Sc  -RL(6,1)*RL(3,2) 

RG(3)=RL(4,2)*RL{3,1)+RL(6,2)*RL(2,1)-RL(4,1)*RL(3,2) 

Sc  -RL(6,1)*RL(2,2) 

RG(4)=RL(4,2)*RL(2,1)+RL(6,2)*RL(1,1)-RL(4, 1)*RL(2,2) 

Sc  -RL(6,1)*RL(1,2)-RL(6,2)*SIGRR+RL(6,1)*SIGTT 

RG(5)=RL(4,2)*RL(1,1)-RL(4,2)*SIGRR-RL(4, 1)*RL(1,2) 

Sc  +RL(4,1)*SIGTT 
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c 

RETURN 

END 

C 

Q  *it*irir1r*ir*irit*1r1tititit*itit'k'k1r*itirit**it1i*itiriciritirie*iciritiritiririt1ririt*iritiririr*iriririt 

SUBROUTINE  NEWTON (RG, RL, SIGRR, SIGRT, A2 , A4 , ATWO, AFOUR, SDIFFl ) 

c 

IMPLICIT  REAL* 8  (A-H,0-Z) 

DIMENSION  RG(6)  ,ATWO{3)  ,AFOUR(3)  .SDIFFIO) 

DIMENSION  RL{10,3) 

C 

TOL=0. 0000000000001 
ITERATION=l 
SDIFF=0 . 0 
A2=0.0 
1=1 

100  FNOT=RG(1)*A2**4.0+RG(2)*A2**3.0+RG(3)*A2**2.0+RG(4)*A2+RG{5) 
STIFF=4.0*RG(1)*A2**3.0+3.0*RG(2) *A2**2.0+2.0*RG{3)*A2+RG(4) 
A2NEW=A2-FNOT/STIFF 
DIFF= ABS ( A2 - A2NEW ) 

IF ( DIFF . GT . TOL ) THEN 

A2=A2NEW 

ITERATIONS ITERATION+ 1 
IF ( ITERATION . GT . 3  0  0 ) THEN 
A2=-1000000000000000000 
A4=-1000000000000000000 
GOTO  200 
ENDIF 
GOTO  100 
ENDIF 
A2=A2NEW 

A4=(SIGRR-(RL(1,1)+RL(2,1)*A2+RL(3,1)*A2**2.0+ 

&  RL(5,1)*A2**3.0) )/(RL(4,l)+RL(6,l)*A2) 

200  ATWO(I)=A2 
AFOUR ( I ) =A4 
C 

SDIFF1{I)=SIGRT-{RL(1,3)+RL(2,3)*A2+RL(3,3) *A2**2.0+RL(4,3)*A4 
&  +RL(5,3)*A2**3.0+RL(6,3) *A2*A4) 

IF ( I . EQ . 1 ) THEN 
1=1+1 
A2=10.0 
GOTO  100 

ELSEIF(I.EQ.2)THEN 
1=1+1 
A2=-10.0 
GOTO  100 

ELSEIF ( I . EQ . 3 ) THEN 

RMAX=10000000000000000 

NUNU=0 

DO  JJ=1,3,1 

IF{ATWO(JJ) .LT.RMAX)THEN 
RMAX=ATWO(JJ) 

NUNU=NUNU+1 

SDIFF=SDIFF1(JJ) 
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A2=ATWO(JJ) 

A4=AFOUR(JJ) 

ENDIF 

ENDDO 

ENDIF 

C 

1  FORMAT{I10,E15.7) 

C 

RETURN 

END 

C 


APPENDIX  H 


FORTRAN  PROGRAM  TO  CALCULATE  STRESSES  FROM 


ANALYTIC  SOLUTIONS 


ooooooooo  o  o  ooo  oo 
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************************************************************** 
PROGRAM  PRINCIPLE_STRESSES 

************************************************************** 


IMPLICIT  NONE 

CHARACTER* 8 0  FILENAME (10). NAME 

INTEGER  NUMFILE,NHARD, J 

INTEGER  TEMP (10) 

REAL*8  R, THETA, DELTHETA 

REAL*8  THETAMAX, ZERO,  PIE,  RT 

REAL*  8  STRESSRTHETA ( 3 ) , DSTRESSRTHETA { 3 ) 

REAL*8  STRESSTEMP(10,2000) ,EIGENCOEF{10,2) ,STRESSK{10, 10) 

REAL*8  STRESS (10, 7, 2000) 

Define  the  data  files  to  be  read  in 


NAME=' eigen. daf 
NUMFILE  =  6 

FILENAME  ( 1 )  = '  f  I_nl0_lainbda5_stress .  dat  * 

FILENAME{2)  =  '  f2_nl0_lainbda5_stress.dat  • 

FILENAME  { 3 )  = '  f  3_nl0_lainbda5_stress . dat  ‘ 

FILENAME  { 4 )  = '  f  4_nl 0_lainbda5_stress .  dat ' 

FILENAME ( 5 ) = ' f 5_nl0_lambda5_stress . dat ’ 

FILENAME  ( 6 )  = '  f  6_nl0_lainbda5_stress .  dat  • 

NHARD=10 

WRITE (*,*) 'Please  enter  the  normalized  radius?' 

read(*, *)R 

R=R*0.002 

THETAMAX=3 .141592654 
DELTHETA=3 . 141592654/90 

PIE=4*ACOS{ {2.0)**{-0.5)) 

RT=360.0/(2.0*PIE) 

ZERO=0. 00000000000000000000000000000000 
Get  angular  stress  functions 

CALL  GETALLSTRESSES  (NUMFILE, FILENAME, STRESS) 

Get  eigen  values  and  coefficients 

CALL  GETEIGENCOE  (NAME, EIGENCOEF, NUMFILE) 

Define  Coefficients 

EIGENCOEF ( 1 , 2 )  =  ( EIGENCOEF (1,2)) *  *EIGENCOEF (1,1) 

EIGENCOEF (3,2) =EIGENCOEF (2 , 2 ) **2 /EIGENCOEF (1,2) 

EIGENCOEF (5,2) =EIGENCOEF ( 2 , 2 ) *  *3 /EIGENCOEF ( 1 , 2 ) *  *2 
EIGENCOEF (6,2) =EIGENCOEF (2,2) *EIGENCOEF (4,2) /EIGENCOEF (1,2) 

Initialize  angle 


C 


THETA=0 
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C  Get  angular  stress  functions  at  theta 
C 

10  CALL  STRESSK_AT_THETA  (THETA, NUMFILE, STRESS, STRESSK) 

C 

C  Get  stress  at  angle 

C 

CALL  STRESS_R_THETA  (NUMFILE, 1, R, EIGENCOEF, STRESSK, 

Sc  STRESSRTHETA,DSTRESSRTHETA) 

C 

C  Write  angle  and  stress  to  output  file 

C 

WRITE(99,1)  THETA, (STRESSRTHETA( J) , J=1 , 3 ) 

C 

C  Increment  theta 

C 

THETA=THETA+DELTHETA 
IF ( THETA. LE.THETAMAX)  GOTO  10 
C 

1  F0RMAT(1X,4(F15.7) ) 

C 

STOP 

END 

C 

SUBROUTINE  STRESS_R_THETA (NUMFILE, BEGIN, R, EIGENCOEF, STRESSK, 

Sc  STRESSRTHETA,  DSTRESSRTHETA) 

C  *************************************************************** 

C 

IMPLICIT  NONE 

INTEGER  I, J, NUMFILE, BEGIN 
REAL*8  R 

REAL*  8  STRESSRTHETA ( 3 ) , DSTRESSRTHETA ( 3 ) 

REAL*8  EIGENCOEF(10,2) ,STRESSK(10,10) 

C 

DO  1=1, 3,1 

STRESSRTHETA ( I ) =  0 
DSTRESSRTHETA ( I ) =  0 
DO  J=BEGIN, NUMFILE,! 

STRESSRTHETA ( I ) =STRESSRTHETA ( 1 ) + 

Sc  EIGENCOEF{J,2)*STRESSK(J,  I)*R**EIGENCOEF(J,l) 

DSTRESSRTHETA ( I ) =DSTRESSRTHETA ( 1 ) + 

Sc  EIGENCOEF  (J,  2  )*EIGENCOEF(J,l)*STRESSK(J,  I)  * 

Sc  R**  (EIGENCOEF (J,l)-1) 

ENDDO 

ENDDO 

C 

RETURN 

END 

C 

^  ************************************************************** 

SUBROUTINE  STRESSK_AT_THETA  (THETA, NUMFILE, STRESS, STRESSK) 

Q  ************************************************************** 

C 

IMPLICIT  NONE 
INTEGER  NUMFILE, I, J,K 


oooo  oo  oo  onoo 
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REAL  *  8  THETA . PERCENT 
REAL*8  STRESSKdO,  10) 

REAL*8  STRESS (10, 7, 2000) 

C 

DO  I=:1,NUMFILE,  1 

DO  J=2, STRESS ( I, 1,1),1 

IF ( (STRESS (1,2, J) .GE. THETA) .AND. 

Sc  (STRESS(I,2,J-1)  .LE. THETA)  )THEN 

PERCENT= (THETA-STRESS ( I , 2 , J-1 ) ) 

Sc  /(STRESS(I,2,J)-STRESS(I,2,  J-1)  ) 

DO  K=l,5,l 

STRESSK (I , K) =STRESS ( I , K+2 , J-1 ) + PERCENT 
Sc  *(STRESS(I,K+2,  J)-STRESS(I,K+2,  J-1)  ) 

ENDDO 
GOTO  1000 
ENDIF 
ENDDO 
1000  ENDDO 
C 

RETURN 
END 

******************'*****'*•'***«'*•*«*«***************'*******««*** 

SUBROUTINE  GETEIGENCOE (NAME, EIGEN, NUMFILE) 

IMPLICIT  NONE 
CHARACTER* 80  NAME 
INTEGER  I, NUMFILE 
REAL*8  EIGEN (10, 2) 

OPEN (UNIT=1 , FILE=NAME, STATUS=  *  OLD • ) 

WRITE (*,*) 'getting  eigenvalues' 

DO  1=1,  NUMFILE,  1 

READd,*)  EIGEN(I,1)  ,EIGEN(I,2) 

ENDDO 

CLOSE  (UNIT=1) 

RETURN 
END 

************************************************************** 

SUBROUTINE  GETALLSTRESSES  (NUMFILE, FILENAME, STRESS) 

*************************************************************** 

IMPLICIT  NONE 

INTEGER  I, J,K,NREC, NUMFILE 

CHARACTER* 80  FILENAME (10) 

REAL* 8  FACTOR 

REAL*8  STRESSTEMPd  0,2000) 

REAL*8  STRESS (10, 7, 2000) 

C 

WRITE (*,*) 'getting  angular  distributions' 
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DO  I=1,NUMFILE, 1 
NREC=0 

CALL  ZERO  ( STRESSTEMP ,10,2000) 

CALL  LAST  (FILENAME { I ), NREC, STRESSTEMP, 6) 

C  CALL  SORT  (NREC, STRESSTEMP, 10, 2000) 

STRESS (1,1,1) =NREC 
DO  .J=2,7,l 

DO  K=1,NREC,1 
FACTOR=1.0 
C 

C  Using  Wang's  convention 

C 

c  IF( (I.GT.l) .AND. (J.GE.3) .AND. (J.LE.6) )THEN 

c  FACTOR=-1.0 

c  ELSE 

c  FACTOR=1.0 

c  ENDIF 

STRESS ( I , J , K ) =FACTOR*STRESSTEMP ( J- 1 , K ) 

ENDDO 

ENDDO 

ENDDO 

C 

RETURN 

END 

C 

Q  W************************************************************* 

SUBROUTINE  LAST (FILENAME, NREC, POTENTIAL, NDIMl ) 

^  n*************************-************************************ 

C 

IMPLICIT  NONE 
CHARACTER* 80  FILENAME 
INTEGER  I, NREC, NDIMl 
REAL*8  POTENTIAL (10, 2000) 

C 

OPEN (UNIT=1 , FILE=FILENAME, STATUS= ' OLD ' ) 

C 

WRITE (*,*)■ reading  in  data' 

NREC  =  1 
READd,  *) 

10  READd,  *,ERR=20,END=20)  (POTENTIAL ( I,  NREC)  ,  1=1,  NDIMl ,  1 ) 

NREC  =  NREC+1 
GOTO  10 

20  NREC  =  NREC-1' 

C 

CLOSE  (UNIT=1) 

C 

RETURN 

END 

C 

Q  *************************************************************** 

SUBROUTINE  SORT (NUM, ARRAY, NDIMl, NDIM2) 

^  ******************•***•*******************«**'****•************ 

C 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  TEMP (10) 


oooonnoo  non 


DIMENSION  ARRAY {NDIM1,NDIM2) 
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Sort  data 

WRITE {*,*) 'sorting' 

DO  I=2,NUM,1 
DO  Jt=2,NUM,  1 

IF ( ARRAY (l.J)  .LT.ARRAYd,  J-1)  )THEN 
DO  K=1,NDIM1,1 

TEMP (K) =ARRAY (K, J-1 ) 

ENDDO 

DO  K=1,NDIM1,1 

ARRAY (K, J-1 ) =ARRAY (K, J) 

ENDDO 

DO  K=1,NDIM1,1 

ARRAY (K,J)=TEMP(K) 

ENDDO 

ENDIF 

ENDDO 

ENDDO 

RETURN 

END 

★★♦*★♦♦★♦★★**★★★*****♦♦★★♦♦♦*★♦*★***♦*★♦♦★*******★*♦★★★★*★★♦*★ 

SUBROUTINE  ZERO (ARRAY, NDIM1,NDIM2) 

********•«*«****«•«««'«********«*•«*********•*•**•««*****••*•** 


IMPLICIT  REAL*8  (A-H,0-Z) 
DIMENSION  ARRAY (NDIM1,NDIM2) 


zero  data 

DO  I=1,NDIM1,1 
DO  J=1,NDIM2,1 
ARRAYd,  J)=0 
ENDDO 
ENDDO 
C 

RETURN 

END 

C 
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